From ebded8e06c53d593946874c468a5d4eff3deb5d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Thu, 16 Oct 2025 10:09:20 +0200 Subject: [PATCH 01/43] Add files via upload MDF taper function for smooth termination of Coulomb integrals (same as used in GULP) --- src/mdftaper.cpp | 23 +++++++++++++++++++++++ src/mdftaper.h | 6 ++++++ 2 files changed, 29 insertions(+) create mode 100644 src/mdftaper.cpp create mode 100644 src/mdftaper.h diff --git a/src/mdftaper.cpp b/src/mdftaper.cpp new file mode 100644 index 00000000000..e940597cb98 --- /dev/null +++ b/src/mdftaper.cpp @@ -0,0 +1,23 @@ +#include + +/* MDF taper function between rmin and rmax as in GULP */ + +void mdftaper(double r, double rmin, double rmax, double &f, double &df) +{ + double x, dx, t1, t2; + + if (r <= rmin) { + f = 1.0; + df = 0.0; + } else if (r >= rmax) { + f = 0.0; + df = 0.0; + } else { + x = (r-rmin)/(rmax-rmin); + dx = 1.0/(rmax-rmin); + t1 = 1.0 - x; + t2 = 1.0 + 3.0*x + 6.0*pow(x,2); + f = pow(t1,3)*t2; + df = (-3.0*pow(t1,2)*t2 + pow(t1,3)*(3.0+12.0*x))*dx; + } +} diff --git a/src/mdftaper.h b/src/mdftaper.h new file mode 100644 index 00000000000..1f2e929ed19 --- /dev/null +++ b/src/mdftaper.h @@ -0,0 +1,6 @@ +#ifndef LMP_TAPER_H +#define LMP_TAPER_H + +void mdftaper(double, double, double, double &, double &); + +#endif From 8154ee26b8064f334471f0567e48b894ac8cb1ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Thu, 16 Oct 2025 10:38:45 +0200 Subject: [PATCH 02/43] Add files via upload Adding smooth taper to Coulomb integrals, allowing calculations using wolf and dsf. --- src/QEQ/fix_qeq_slater.cpp | 106 ++++++++++++++++++++++++++++++++++--- src/QEQ/fix_qeq_slater.h | 3 ++ 2 files changed, 102 insertions(+), 7 deletions(-) diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index e80bad04b97..67afc14e663 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -23,6 +23,7 @@ #include "force.h" #include "kspace.h" #include "math_const.h" +#include "mdftaper.h" #include "neigh_list.h" #include "neighbor.h" #include "pair.h" @@ -40,6 +41,8 @@ using namespace FixConst; FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg, arg) { alpha = 0.20; + vtype = 0; + drtap = 0.0; // optional arg int iarg = 8; @@ -48,12 +51,31 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/slater alpha", error); alpha = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; + if (iarg < narg) { + if (strcmp(arg[iarg],"wolf") == 0) { // type of potential (0=erfc(r)/r, 1=Wolf, 2=Fennell & Gezelter) + vtype = 1; + iarg++; + if (iarg < narg) { + drtap = utils::numeric(FLERR, arg[iarg], false, lmp); + iarg++; + } + } else if (strcmp(arg[iarg],"dsf") == 0) { + vtype = 2; + iarg++; + if (iarg < narg) { + drtap = utils::numeric(FLERR, arg[iarg], false, lmp); + iarg++; + } + } + } } else if (strcmp(arg[iarg], "warn") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/slater warn", error); maxwarn = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else + } else { error->all(FLERR, "Unknown fix qeq/slater keyword: {}", arg[iarg]); + continue; + } } if (streitz_flag) extract_streitz(); @@ -189,7 +211,12 @@ void FixQEqSlater::compute_H() r = sqrt(rsq); H.jlist[m_fill] = j; - H.val[m_fill] = calculate_H(zei, zej, zj, r, zjtmp); + if (vtype == 0) + H.val[m_fill] = calculate_H(zei, zej, zj, r, zjtmp); + else if (vtype == 1) + H.val[m_fill] = calculate_H_wolf(zei, zej, zj, r, zjtmp); + else if (vtype == 2) + H.val[m_fill] = calculate_H_dsf(zei, zej, zj, r, zjtmp); m_fill++; } H.numnbrs[i] = m_fill - H.firstnbr[i]; @@ -207,6 +234,7 @@ double FixQEqSlater::calculate_H(double zei, double zej, double zj, double r, double &zjtmp) { double rinv = 1.0/r; + double rc = cutoff; double exp2zir = exp(-2.0*zei*r); double zei2 = zei*zei; @@ -228,6 +256,7 @@ double FixQEqSlater::calculate_H(double zei, double zej, double zj, double etmp1, etmp2; double e1, e2, e3, e4; double ci_jfi, ci_fifj; + double ftap, dftap; e1 = e2 = e3 = e4 = 0.0; etmp1 = etmp2 = 0.0; @@ -246,8 +275,10 @@ double FixQEqSlater::calculate_H(double zei, double zej, double zj, ci_fifj = -exp2zir*(e1+e3/r) - exp2zjr*(e2+e4/r); } - etmp1 = 1.00 * (ci_jfi - ci_fifj); - etmp2 = 0.50 * (ci_fifj + erfcr*rinv); + mdftaper(r, rc-drtap, rc, ftap, dftap); + + etmp1 = 1.00 * (ci_jfi - ci_fifj) * ftap; + etmp2 = 0.50 * (ci_fifj + erfcr*rinv) * ftap; zjtmp += qqrd2e * zj * etmp1; return qqrd2e * etmp2; @@ -284,6 +315,7 @@ double FixQEqSlater::calculate_H_wolf(double zei, double zej, double zj, double eshift, fshift, ci_jfi, ci_fifj; double etmp1, etmp2, etmp3; + double ftap, dftap; double a = alpha; double erfcr = erfc(a*r); @@ -316,13 +348,74 @@ double FixQEqSlater::calculate_H_wolf(double zei, double zej, double zj, - eshift - (r-rc)*fshift; } + mdftaper(r, rc-drtap, rc, ftap, dftap); + etmp1 = erfcr/r - erfcrc/rc; - etmp2 = 1.00 * (ci_jfi - ci_fifj); - etmp3 = 0.50 * (etmp1 + ci_fifj); + etmp2 = 1.00 * (ci_jfi - ci_fifj) * ftap; + etmp3 = 0.50 * (etmp1 + ci_fifj*ftap); zjtmp += qqrd2e * zj * etmp2; return qqrd2e * etmp3; +} + +/* ---------------------------------------------------------------------- */ + +double FixQEqSlater::calculate_H_dsf(double zei, double zej, double zj, + double r, double &zjtmp) +{ + double rinv = 1.0/r; + + double exp2zir = exp(-2.0*zei*r); + double zei2 = zei*zei; + double zei4 = zei2*zei2; + double zei6 = zei2*zei4; + + double exp2zjr = exp(-2.0*zej*r); + double zej2 = zej*zej; + double zej4 = zej2*zej2; + double zej6 = zej2*zej4; + + double sm1 = 11.0/8.0; + double sm2 = 3.00/4.0; + double sm3 = 1.00/6.0; + + double erfcr = erfc(alpha*r); + double qqrd2e = force->qqrd2e; + double etmp1, etmp2; + double e1, e2, e3, e4; + double ci_jfi, ci_fifj; + + double rc = cutoff; + double erfcrc = erfc(alpha*rc); + double edsf; + double ftap, dftap; + + e1 = e2 = e3 = e4 = 0.0; + etmp1 = etmp2 = edsf = 0.0; + + ci_jfi = -zei*exp2zir - rinv*exp2zir; + + if (zei == zej) { + ci_fifj = -exp2zir*(rinv + zei*(sm1 + sm2*zei*r + sm3*zei2*r*r)); + } else { + e1 = zei*zej4/((zei+zej)*(zei+zej)*(zei-zej)*(zei-zej)); + e2 = zej*zei4/((zei+zej)*(zei+zej)*(zej-zei)*(zej-zei)); + e3 = (3.0*zei2*zej4-zej6) / + ((zei+zej)*(zei+zej)*(zei+zej)*(zei-zej)*(zei-zej)*(zei-zej)); + e4 = (3.0*zej2*zei4-zei6) / + ((zei+zej)*(zei+zej)*(zei+zej)*(zej-zei)*(zej-zei)*(zej-zei)); + ci_fifj = -exp2zir*(e1+e3/r) - exp2zjr*(e2+e4/r); + } + + mdftaper(r, rc-drtap, rc, ftap, dftap); + + edsf = erfcr/r - erfcrc/rc + (erfcrc/(rc*rc)+2.0*alpha/MY_PIS*exp(-alpha*alpha*rc*rc)/rc)*(r-rc); + etmp1 = 1.00 * (ci_jfi - ci_fifj) * ftap; + etmp2 = 0.50 * (ci_fifj*ftap + edsf); + + zjtmp += qqrd2e * zj * etmp1; + return qqrd2e * etmp2; } /* ---------------------------------------------------------------------- */ @@ -356,5 +449,4 @@ void FixQEqSlater::sparse_matvec(sparse_matrix *A, double *x, double *b) } } } - } diff --git a/src/QEQ/fix_qeq_slater.h b/src/QEQ/fix_qeq_slater.h index 69b31c8066f..94a6bef0cfc 100644 --- a/src/QEQ/fix_qeq_slater.h +++ b/src/QEQ/fix_qeq_slater.h @@ -37,10 +37,13 @@ class FixQEqSlater : public FixQEq { void compute_H(); double calculate_H(double, double, double, double, double &); double calculate_H_wolf(double, double, double, double, double &); + double calculate_H_dsf(double, double, double, double, double &); void extract_streitz(); class PairCoulStreitz *streitz; double alpha; + double drtap; + int vtype; }; } // namespace LAMMPS_NS #endif From 95a88634db5f6445490fdf29a85b6bb7906aeff4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Thu, 16 Oct 2025 10:39:41 +0200 Subject: [PATCH 03/43] Add files via upload Adding smooth taper to Coulomb integrals. --- src/KSPACE/pair_coul_streitz.cpp | 86 ++++++++++++++++++++++++++++---- src/KSPACE/pair_coul_streitz.h | 5 ++ 2 files changed, 81 insertions(+), 10 deletions(-) diff --git a/src/KSPACE/pair_coul_streitz.cpp b/src/KSPACE/pair_coul_streitz.cpp index 055d98309bf..86a0291135c 100644 --- a/src/KSPACE/pair_coul_streitz.cpp +++ b/src/KSPACE/pair_coul_streitz.cpp @@ -25,6 +25,7 @@ #include "info.h" #include "kspace.h" #include "math_const.h" +#include "mdftaper.h" #include "memory.h" #include "neigh_list.h" #include "neighbor.h" @@ -47,6 +48,7 @@ PairCoulStreitz::PairCoulStreitz(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; nmax = 0; + drtap = 0.0; params = nullptr; } @@ -102,14 +104,23 @@ void PairCoulStreitz::settings(int narg, char **arg) cut_coul = utils::numeric(FLERR,arg[0],false,lmp); if (strcmp(arg[1],"wolf") == 0) { + if (narg < 3) utils::missing_cmd_args(FLERR, "pair_style coul/streitz wolf", error); + kspacetype = 1; + g_wolf = utils::numeric(FLERR,arg[2],false,lmp); + if (narg > 3) + drtap = utils::numeric(FLERR,arg[3],false,lmp); + } else if (strcmp(arg[1],"dsf") == 0){ + if (narg < 3) utils::missing_cmd_args(FLERR, "pair_style coul/streitz dsf", error); kspacetype = 1; + dsfflag = 1; g_wolf = utils::numeric(FLERR,arg[2],false,lmp); + if (narg > 3) + drtap = utils::numeric(FLERR,arg[3],false,lmp); } else if (strcmp(arg[1],"ewald") == 0) { ewaldflag = pppmflag = 1; kspacetype = 2; - } else { - error->all(FLERR,"Illegal pair_style command"); - } + } else + error->all(FLERR, "Unknown pair_style keyword: {}", arg[1]); } /* ---------------------------------------------------------------------- @@ -358,10 +369,14 @@ void PairCoulStreitz::compute(int eflag, int vflag) coulomb_integral_wolf(zei, zej, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj); - // Wolf Sum + // Wolf-Fennell Sum - wolf_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, - ecoul, forcecoul); + if (dsfflag == 1) + fennell_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, + ecoul, forcecoul); + else + wolf_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, + ecoul, forcecoul); // Forces @@ -558,20 +573,71 @@ void PairCoulStreitz::wolf_sum(double qi, double qj, double zj, double r, double etmp1, etmp2, etmp3; double ftmp1, ftmp2, ftmp3; + double ftap, dftap; + etmp = etmp1 = etmp2 = etmp3 = 0.0; ftmp = ftmp1 = ftmp2 = ftmp3 = 0.0; + mdftaper(r, rc-drtap, rc, ftap, dftap); + etmp1 = erfcr/r - erfcrc/rc; - etmp2 = qi * zj * (ci_jfi - ci_fifj); - etmp3 = qi * qj * 0.50 * (etmp1 + ci_fifj); + etmp2 = qi * zj * (ci_jfi - ci_fifj) * ftap; + etmp3 = qi * qj * 0.50 * (etmp1 + ci_fifj*ftap); ftmp1 = -erfcr/r/r - 2.0*a/MY_PIS*derfcr/r - dwoself; - ftmp2 = qi * zj * (dci_jfi - dci_fifj); - ftmp3 = qi * qj * 0.50 * (ftmp1 + dci_fifj); + ftmp2 = qi * zj * (dci_jfi*ftap + ci_jfi*dftap - dci_fifj*ftap - ci_fifj*dftap); + ftmp3 = qi * qj * 0.50 * (ftmp1 + dci_fifj*ftap + ci_fifj*dftap); etmp = qqrd2e * (etmp2 + etmp3); ftmp = qqrd2e * (ftmp2 + ftmp3); +} + +/* ---------------------------------------------------------------------- + zi*zj terms commented out, because they are contained in the potential (pair, EAM, Tersoff, ...) +*/ + +void PairCoulStreitz::fennell_sum(double qi, double qj, double zj, double r, + double ci_jfi, double dci_jfi, double ci_fifj, + double dci_fifj, double &etmp, double &ftmp) +{ + double a = g_wolf; + double rc = cut_coul; + double qqrd2e = force->qqrd2e; + + double r2 = r*r, rc2 = rc*rc, a2 = a*a; + double erfcr = erfc(a*r); + double derfcr = exp(-a2*r2); + double derfcrc = exp(-a2*rc2); + double erfcrc = erfc(a*rc); + + double etmp0, /*etmp1,*/ etmp2, etmp3; + double ftmp0, /*ftmp1,*/ ftmp2, ftmp3; + double con = 2.0*a/MY_PIS; + double ftap, dftap; + + if (r >= rc) { + etmp = 0.0; + ftmp = 0.0; + return; + } + + mdftaper(r, rc-drtap, rc, ftap, dftap); + + etmp0 = erfcr/r - erfcrc/rc + (erfcrc/rc2 + con*derfcrc/rc)*(r-rc); +// etmp1 = zi * zj * 0.5 * (ci_fifj*ftap - ci_ifj*ftap - ci_jfi*ftap); + etmp2 = qi * zj * (ci_jfi - ci_fifj) * ftap; + etmp3 = qi * qj * 0.5 * (etmp0 + ci_fifj*ftap); + + ftmp0 = -erfcr/r2 - con*derfcr/r + (erfcrc/rc2 + con*derfcrc/rc); + ftmp0 = ftmp0 - dwoself; +// ftmp1 = zi * zj * 0.5 * (dci_fifj*ftap + ci_fifj*dftap - dci_ifj*ftap - ci_ifj*dftap - +// dci_jfi*ftap - ci_jfi*dftap); + ftmp2 = qi * zj * (dci_jfi*ftap + ci_jfi*dftap - dci_fifj*ftap - ci_fifj*dftap); + ftmp3 = qi * qj * 0.5 * (ftmp0 + dci_fifj*ftap + ci_fifj*dftap); + + etmp = qqrd2e * (/*etmp1 +*/ etmp2 + etmp3); + ftmp = qqrd2e * (/*ftmp1 +*/ ftmp2 + ftmp3); } /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_coul_streitz.h b/src/KSPACE/pair_coul_streitz.h index 059ec3f6ae2..2fb8c8041eb 100644 --- a/src/KSPACE/pair_coul_streitz.h +++ b/src/KSPACE/pair_coul_streitz.h @@ -56,6 +56,10 @@ class PairCoulStreitz : public Pair { // Wolf double g_wolf, woself, dwoself; + // Fennell-Gezelter + double drtap; + int dsfflag; + // Ewald double g_ewald; @@ -68,6 +72,7 @@ class PairCoulStreitz : public Pair { double self(Param *, double); void coulomb_integral_wolf(double, double, double, double &, double &, double &, double &); void wolf_sum(double, double, double, double, double, double, double, double, double &, double &); + void fennell_sum(double, double, double, double, double, double, double, double, double &, double &); void coulomb_integral_ewald(double, double, double, double &, double &, double &, double &); void ewald_sum(double, double, double, double, double, double, double, double, double &, double &, double); From 2b5144cf7fae375f2fdabbf0d789efd020caefdc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Fri, 17 Oct 2025 09:07:00 +0200 Subject: [PATCH 04/43] Add files via upload --- src/QEQ/fix_qeq_slater.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index 67afc14e663..fb50671143a 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -23,7 +23,7 @@ #include "force.h" #include "kspace.h" #include "math_const.h" -#include "mdftaper.h" +#include "math_special.h" #include "neigh_list.h" #include "neighbor.h" #include "pair.h" @@ -35,6 +35,8 @@ using namespace LAMMPS_NS; using namespace MathConst; using namespace FixConst; +using MathSpecial::mdftaper; + /* ---------------------------------------------------------------------- */ From 704f7f822c718fc0b528c8a29e076d94c8643282 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Fri, 17 Oct 2025 09:08:01 +0200 Subject: [PATCH 05/43] Add files via upload --- src/KSPACE/pair_coul_streitz.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/KSPACE/pair_coul_streitz.cpp b/src/KSPACE/pair_coul_streitz.cpp index 86a0291135c..e30ca49756a 100644 --- a/src/KSPACE/pair_coul_streitz.cpp +++ b/src/KSPACE/pair_coul_streitz.cpp @@ -25,7 +25,7 @@ #include "info.h" #include "kspace.h" #include "math_const.h" -#include "mdftaper.h" +#include "math_special.h" #include "memory.h" #include "neigh_list.h" #include "neighbor.h" @@ -36,6 +36,7 @@ using namespace LAMMPS_NS; using namespace MathConst; +using MathSpecial::mdftaper; static constexpr int DELTA = 4; static constexpr int MAXNEIGH = 24; From f08c3d3c5320aedf84b562e0e385d0cfa792ae31 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Fri, 17 Oct 2025 09:10:00 +0200 Subject: [PATCH 06/43] Add files via upload Added function mdftaper() --- src/math_special.cpp | 20 ++++++++++++++++++++ src/math_special.h | 5 +++++ 2 files changed, 25 insertions(+) diff --git a/src/math_special.cpp b/src/math_special.cpp index e7cd4ef75d8..c23cb23228c 100644 --- a/src/math_special.cpp +++ b/src/math_special.cpp @@ -745,3 +745,23 @@ double MathSpecial::fm_exp(double x) #endif } +void MathSpecial::mdftaper(double r, double rmin, double rmax, double &f, double &df) +{ + double x, dx, t1, t2, cube_t1; + + if (r <= rmin) { + f = 1.0; + df = 0.0; + } else if (r >= rmax) { + f = 0.0; + df = 0.0; + } else { + x = (r-rmin)/(rmax-rmin); + dx = 1.0/(rmax-rmin); + t1 = 1.0 - x; + t2 = 1.0 + 3.0*x + 6.0*square(x); + cube_t1 = cube(t1); + f = cube_t1*t2; + df = (-3.0*square(t1)*t2 + cube_t1*(3.0+12.0*x))*dx; + } +} diff --git a/src/math_special.h b/src/math_special.h index b190a4b622d..11eda2569a7 100644 --- a/src/math_special.h +++ b/src/math_special.h @@ -60,6 +60,11 @@ namespace MathSpecial { extern double fm_exp(double x); + /* MDF taper function defined in Mei et al. (Phys. Rev. B 43:4653, 1991). It is used + * to smoothly terminate functions between rmin and rmax. */ + + extern void mdftaper(double r, double rmin, double rmax, double &f, double &df); + // support function for scaled error function complement extern double erfcx_y100(const double y100); From 1e114dcfb48b1c313b080a022cf5de0b3f20bbfb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Fri, 17 Oct 2025 09:11:17 +0200 Subject: [PATCH 07/43] Delete src/mdftaper.cpp mdftaper() moved into math_special --- src/mdftaper.cpp | 23 ----------------------- 1 file changed, 23 deletions(-) delete mode 100644 src/mdftaper.cpp diff --git a/src/mdftaper.cpp b/src/mdftaper.cpp deleted file mode 100644 index e940597cb98..00000000000 --- a/src/mdftaper.cpp +++ /dev/null @@ -1,23 +0,0 @@ -#include - -/* MDF taper function between rmin and rmax as in GULP */ - -void mdftaper(double r, double rmin, double rmax, double &f, double &df) -{ - double x, dx, t1, t2; - - if (r <= rmin) { - f = 1.0; - df = 0.0; - } else if (r >= rmax) { - f = 0.0; - df = 0.0; - } else { - x = (r-rmin)/(rmax-rmin); - dx = 1.0/(rmax-rmin); - t1 = 1.0 - x; - t2 = 1.0 + 3.0*x + 6.0*pow(x,2); - f = pow(t1,3)*t2; - df = (-3.0*pow(t1,2)*t2 + pow(t1,3)*(3.0+12.0*x))*dx; - } -} From 0b2e03d3de1edf7955646a6fb5678875aefe86e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Fri, 17 Oct 2025 09:11:36 +0200 Subject: [PATCH 08/43] Delete src/mdftaper.h mdftaper() moved into math_special --- src/mdftaper.h | 6 ------ 1 file changed, 6 deletions(-) delete mode 100644 src/mdftaper.h diff --git a/src/mdftaper.h b/src/mdftaper.h deleted file mode 100644 index 1f2e929ed19..00000000000 --- a/src/mdftaper.h +++ /dev/null @@ -1,6 +0,0 @@ -#ifndef LMP_TAPER_H -#define LMP_TAPER_H - -void mdftaper(double, double, double, double &, double &); - -#endif From 365fbe4ab678df2445e998815c5e1038e7cab6b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Wed, 22 Oct 2025 13:48:38 +0200 Subject: [PATCH 09/43] Add files via upload Removed trailing spaces, converted tabs to spaces --- src/QEQ/fix_qeq_slater.cpp | 40 +++++++++++++++++++------------------- src/QEQ/fix_qeq_slater.h | 4 ++-- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index fb50671143a..be0e2346ab3 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -54,21 +54,21 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg alpha = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; if (iarg < narg) { - if (strcmp(arg[iarg],"wolf") == 0) { // type of potential (0=erfc(r)/r, 1=Wolf, 2=Fennell & Gezelter) - vtype = 1; - iarg++; - if (iarg < narg) { - drtap = utils::numeric(FLERR, arg[iarg], false, lmp); - iarg++; - } - } else if (strcmp(arg[iarg],"dsf") == 0) { - vtype = 2; - iarg++; - if (iarg < narg) { - drtap = utils::numeric(FLERR, arg[iarg], false, lmp); - iarg++; - } - } + if (strcmp(arg[iarg],"wolf") == 0) { // type of potential (0=erfc(r)/r, 1=Wolf, 2=Fennell & Gezelter) + vtype = 1; + iarg++; + if (iarg < narg) { + drtap = utils::numeric(FLERR, arg[iarg], false, lmp); + iarg++; + } + } else if (strcmp(arg[iarg],"dsf") == 0) { + vtype = 2; + iarg++; + if (iarg < narg) { + drtap = utils::numeric(FLERR, arg[iarg], false, lmp); + iarg++; + } + } } } else if (strcmp(arg[iarg], "warn") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix qeq/slater warn", error); @@ -214,11 +214,11 @@ void FixQEqSlater::compute_H() r = sqrt(rsq); H.jlist[m_fill] = j; if (vtype == 0) - H.val[m_fill] = calculate_H(zei, zej, zj, r, zjtmp); + H.val[m_fill] = calculate_H(zei, zej, zj, r, zjtmp); else if (vtype == 1) - H.val[m_fill] = calculate_H_wolf(zei, zej, zj, r, zjtmp); + H.val[m_fill] = calculate_H_wolf(zei, zej, zj, r, zjtmp); else if (vtype == 2) - H.val[m_fill] = calculate_H_dsf(zei, zej, zj, r, zjtmp); + H.val[m_fill] = calculate_H_dsf(zei, zej, zj, r, zjtmp); m_fill++; } H.numnbrs[i] = m_fill - H.firstnbr[i]; @@ -236,7 +236,7 @@ double FixQEqSlater::calculate_H(double zei, double zej, double zj, double r, double &zjtmp) { double rinv = 1.0/r; - double rc = cutoff; + double rc = cutoff; double exp2zir = exp(-2.0*zei*r); double zei2 = zei*zei; @@ -317,7 +317,7 @@ double FixQEqSlater::calculate_H_wolf(double zei, double zej, double zj, double eshift, fshift, ci_jfi, ci_fifj; double etmp1, etmp2, etmp3; - double ftap, dftap; + double ftap, dftap; double a = alpha; double erfcr = erfc(a*r); diff --git a/src/QEQ/fix_qeq_slater.h b/src/QEQ/fix_qeq_slater.h index 94a6bef0cfc..f84fa8e417c 100644 --- a/src/QEQ/fix_qeq_slater.h +++ b/src/QEQ/fix_qeq_slater.h @@ -37,13 +37,13 @@ class FixQEqSlater : public FixQEq { void compute_H(); double calculate_H(double, double, double, double, double &); double calculate_H_wolf(double, double, double, double, double &); - double calculate_H_dsf(double, double, double, double, double &); + double calculate_H_dsf(double, double, double, double, double &); void extract_streitz(); class PairCoulStreitz *streitz; double alpha; double drtap; - int vtype; + int vtype; }; } // namespace LAMMPS_NS #endif From b46e4147ad6c3d6ac73b5ce691c4369ae4de9050 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Wed, 22 Oct 2025 13:49:25 +0200 Subject: [PATCH 10/43] Add files via upload Removed trailing spaces, converted tabs to spaces --- src/KSPACE/pair_coul_streitz.cpp | 22 +++++++++++----------- src/KSPACE/pair_coul_streitz.h | 2 +- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/KSPACE/pair_coul_streitz.cpp b/src/KSPACE/pair_coul_streitz.cpp index e30ca49756a..9206c41796c 100644 --- a/src/KSPACE/pair_coul_streitz.cpp +++ b/src/KSPACE/pair_coul_streitz.cpp @@ -106,16 +106,16 @@ void PairCoulStreitz::settings(int narg, char **arg) if (strcmp(arg[1],"wolf") == 0) { if (narg < 3) utils::missing_cmd_args(FLERR, "pair_style coul/streitz wolf", error); - kspacetype = 1; + kspacetype = 1; g_wolf = utils::numeric(FLERR,arg[2],false,lmp); if (narg > 3) drtap = utils::numeric(FLERR,arg[3],false,lmp); } else if (strcmp(arg[1],"dsf") == 0){ if (narg < 3) utils::missing_cmd_args(FLERR, "pair_style coul/streitz dsf", error); kspacetype = 1; - dsfflag = 1; + dsfflag = 1; g_wolf = utils::numeric(FLERR,arg[2],false,lmp); - if (narg > 3) + if (narg > 3) drtap = utils::numeric(FLERR,arg[3],false,lmp); } else if (strcmp(arg[1],"ewald") == 0) { ewaldflag = pppmflag = 1; @@ -373,11 +373,11 @@ void PairCoulStreitz::compute(int eflag, int vflag) // Wolf-Fennell Sum if (dsfflag == 1) - fennell_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, - ecoul, forcecoul); - else - wolf_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, - ecoul, forcecoul); + fennell_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, + ecoul, forcecoul); + else + wolf_sum(qi, qj, zj, r, ci_jfi, dci_jfi, ci_fifj, dci_fifj, + ecoul, forcecoul); // Forces @@ -580,7 +580,7 @@ void PairCoulStreitz::wolf_sum(double qi, double qj, double zj, double r, ftmp = ftmp1 = ftmp2 = ftmp3 = 0.0; mdftaper(r, rc-drtap, rc, ftap, dftap); - + etmp1 = erfcr/r - erfcrc/rc; etmp2 = qi * zj * (ci_jfi - ci_fifj) * ftap; etmp3 = qi * qj * 0.50 * (etmp1 + ci_fifj*ftap); @@ -622,14 +622,14 @@ void PairCoulStreitz::fennell_sum(double qi, double qj, double zj, double r, ftmp = 0.0; return; } - + mdftaper(r, rc-drtap, rc, ftap, dftap); etmp0 = erfcr/r - erfcrc/rc + (erfcrc/rc2 + con*derfcrc/rc)*(r-rc); // etmp1 = zi * zj * 0.5 * (ci_fifj*ftap - ci_ifj*ftap - ci_jfi*ftap); etmp2 = qi * zj * (ci_jfi - ci_fifj) * ftap; etmp3 = qi * qj * 0.5 * (etmp0 + ci_fifj*ftap); - + ftmp0 = -erfcr/r2 - con*derfcr/r + (erfcrc/rc2 + con*derfcrc/rc); ftmp0 = ftmp0 - dwoself; // ftmp1 = zi * zj * 0.5 * (dci_fifj*ftap + ci_fifj*dftap - dci_ifj*ftap - ci_ifj*dftap - diff --git a/src/KSPACE/pair_coul_streitz.h b/src/KSPACE/pair_coul_streitz.h index 2fb8c8041eb..7e386cecbf1 100644 --- a/src/KSPACE/pair_coul_streitz.h +++ b/src/KSPACE/pair_coul_streitz.h @@ -72,7 +72,7 @@ class PairCoulStreitz : public Pair { double self(Param *, double); void coulomb_integral_wolf(double, double, double, double &, double &, double &, double &); void wolf_sum(double, double, double, double, double, double, double, double, double &, double &); - void fennell_sum(double, double, double, double, double, double, double, double, double &, double &); + void fennell_sum(double, double, double, double, double, double, double, double, double &, double &); void coulomb_integral_ewald(double, double, double, double &, double &, double &, double &); void ewald_sum(double, double, double, double, double, double, double, double, double &, double &, double); From cdf136d346482179917df412e175b078f0fcd0de Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 14 Nov 2025 10:09:22 -0700 Subject: [PATCH 11/43] Updating dw and arg index bugs --- src/RHEO/compute_rheo_vshift.cpp | 3 ++- src/RHEO/fix_rheo.cpp | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 033d34f899f..135950c9026 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -318,7 +318,7 @@ void ComputeRHEOVShift::correct_type_interface() { int i, j, a, ii, jj, jnum, itype, jtype; int fluidi, fluidj; - double xtmp, ytmp, ztmp, rsq, r, w; + double xtmp, ytmp, ztmp, rsq, r, w, wp; double imass, jmass, voli, volj, rhoi, rhoj; double dx[3]; int dim = domain->dimension; @@ -470,6 +470,7 @@ void ComputeRHEOVShift::correct_type_interface() volj = jmass / rhoj; w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r); + wp = compute_kernel->calc_wp(i, j, dx[0], dx[1], dx[2], r); dWij = compute_kernel->dWij; dWji = compute_kernel->dWji; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 88dd90cdc15..37d58f656bf 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -159,7 +159,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg], "rho/sum") == 0) { rhosum_flag = 1; while (iarg < narg) { // optional sub-arguments - if (strcmp(arg[iarg], "self/mass") == 0) { + if (strcmp(arg[iarg + 1], "self/mass") == 0) { rhosum_self_mass_flag = 1; } else { break; From b6481d5f83be63a20b33b46a17d983f2073302b9 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 14 Nov 2025 12:07:19 -0700 Subject: [PATCH 12/43] Typos in comments --- src/RHEO/compute_rheo_kernel.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 4778eb7b770..c7db53585a8 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -714,7 +714,7 @@ void ComputeRHEOKernel::compute_peratom() if (coordination[i] < zmin) continue; // Use LAPACK to get Minv, use Cholesky decomposition since the - // polynomials are independent, M is symmetrix & positive-definite + // polynomials are independent, M is symmetric & positive-definite const char uplo = 'U'; dpotrf_(&uplo, &Mdim, M, &Mdim, &lapack_error); @@ -742,7 +742,7 @@ void ComputeRHEOKernel::compute_peratom() } // Correction coefficients are columns of M^-1 multiplied by an appropriate coefficient - // Solve the linear system several times to get coefficientns + // Solve the linear system several times to get coefficients // M: 1 x y (z) x^2 y^2 (z^2) xy (xz) (yz) // ---------------------------------------------------------- // 0 1 2 3 4 5 || 2D indexing @@ -764,7 +764,7 @@ void ComputeRHEOKernel::compute_peratom() for (b = 0; b < dim; b++) { //First derivatives C[i][1 + b][a] = -M[a * Mdim + b + 1] * cutinv; - // columns 1-2 (2D) or 1-3 (3D) + // columns 1-2 (2D) or 1-3 (3D) //Second derivatives if (kernel_style == RK2) C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * cutsqinv; From 0ad805f95c3f16139cb283ac6c27c663a58931ff Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 14 Nov 2025 12:08:20 -0700 Subject: [PATCH 13/43] Correcting that wmin criterion is actually a weighted distance --- src/RHEO/compute_rheo_vshift.cpp | 24 ++++++++++++------------ src/RHEO/compute_rheo_vshift.h | 4 ++-- src/RHEO/fix_rheo.cpp | 2 +- src/RHEO/fix_rheo.h | 2 +- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 135950c9026..7aad827edfc 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -42,7 +42,7 @@ using namespace MathExtra; /* ---------------------------------------------------------------------- */ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), wsame(nullptr), + Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), rsame(nullptr), ct(nullptr), cgradt(nullptr), shift_type(nullptr), list(nullptr), compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) { @@ -84,7 +84,7 @@ void ComputeRHEOVShift::init() cross_type_flag = fix_rheo->shift_cross_type_flag; if (cross_type_flag) { scale = fix_rheo->shift_scale; - wmin = fix_rheo->shift_wmin; + rmin = fix_rheo->shift_rmin; cmin = fix_rheo->shift_cmin; comm_forward = 1; comm_reverse = 4; @@ -137,7 +137,7 @@ void ComputeRHEOVShift::compute_peratom() if (cross_type_flag) { memory->grow(ct, atom->nmax, "rheo:ct"); memory->grow(cgradt, atom->nmax, 3, "rheo:cgradt"); - memory->grow(wsame, atom->nmax, "rheo:wsame"); + memory->grow(rsame, atom->nmax, "rheo:rsame"); } nmax_store = atom->nmax; @@ -343,7 +343,7 @@ void ComputeRHEOVShift::correct_type_interface() size_t nbytes = nmax_store * sizeof(double); memset(&ct[0], 0, nbytes); - memset(&wsame[0], 0, nbytes); + memset(&rsame[0], 0, nbytes); memset(&cgradt[0][0], 0, 3 * nbytes); double ctmp, *dWij, *dWji; @@ -470,7 +470,7 @@ void ComputeRHEOVShift::correct_type_interface() volj = jmass / rhoj; w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r); - wp = compute_kernel->calc_wp(i, j, dx[0], dx[1], dx[2], r); + wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2], r); dWij = compute_kernel->dWij; dWji = compute_kernel->dWji; @@ -484,9 +484,9 @@ void ComputeRHEOVShift::correct_type_interface() } if (itype == jtype) { - wsame[i] += w * r; + rsame[i] += w * r; if (newton_pair || j < nlocal) - wsame[j] += w * r; + rsame[j] += w * r; } } } @@ -504,7 +504,7 @@ void ComputeRHEOVShift::correct_type_interface() for (i = 0; i < nlocal; i++) { // If isolated, just don't shift - if (wsame[i] < wmin) { + if (rsame[i] < rmin) { for (a = 0; a < dim; a++) vshift[i][a] = 0.0; continue; @@ -543,7 +543,7 @@ int ComputeRHEOVShift::pack_forward_comm(int n, int *list, double *buf, int /*pb for (i = 0; i < n; i++) { j = list[i]; - buf[m++] = wsame[j]; + buf[m++] = rsame[j]; } return m; @@ -558,7 +558,7 @@ void ComputeRHEOVShift::unpack_forward_comm(int n, int first, double *buf) m = 0; last = first + n; for (i = first; i < last; i++) - wsame[i] = buf[m++]; + rsame[i] = buf[m++]; } /* ---------------------------------------------------------------------- */ @@ -582,7 +582,7 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) for (i = first; i < last; i++) { for (a = 0; a < 3; a++) buf[m++] = cgradt[i][a]; - buf[m++] = wsame[i]; + buf[m++] = rsame[i]; } } return m; @@ -612,7 +612,7 @@ void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf) j = list[i]; for (a = 0; a < 3; a++) cgradt[j][a] += buf[m++]; - wsame[j] += buf[m++]; + rsame[j] += buf[m++]; } } } diff --git a/src/RHEO/compute_rheo_vshift.h b/src/RHEO/compute_rheo_vshift.h index b76d48d5e3d..2913bdd61e8 100644 --- a/src/RHEO/compute_rheo_vshift.h +++ b/src/RHEO/compute_rheo_vshift.h @@ -45,10 +45,10 @@ class ComputeRHEOVShift : public Compute { private: int nmax_store, comm_stage; double dtv, cut, cutsq, cutthird; - double scale, wmin, cmin; + double scale, rmin, cmin; int surface_flag, interface_flag, cross_type_flag; double *rho0; - double *wsame, *ct, **cgradt; + double *rsame, *ct, **cgradt; int *shift_type; class NeighList *list; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index 37d58f656bf..a9640bafe94 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -125,7 +125,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : shift_cross_type_flag = 1; shift_scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - shift_wmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + shift_rmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp); iarg += 3; } else if (strcmp(arg[iarg + 1], "exclude/type") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index aa80bc938e2..a0243c92770 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -57,7 +57,7 @@ class FixRHEO : public Fix { int rhosum_self_mass_flag; int *shift_type; int shift_cross_type_flag; - double shift_scale, shift_wmin, shift_cmin; + double shift_scale, shift_rmin, shift_cmin; // Accessory fixes/computes int viscosity_fix_defined; From 84328a1420de798969e806af307179aa8509ffda Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 14 Nov 2025 14:09:12 -0700 Subject: [PATCH 14/43] Updating wmin -> rmin in doc files --- doc/src/fix_rheo.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/src/fix_rheo.rst b/doc/src/fix_rheo.rst index ef18d19f82e..ba4452fd46b 100644 --- a/doc/src/fix_rheo.rst +++ b/doc/src/fix_rheo.rst @@ -34,10 +34,10 @@ Syntax optional args = *exclude/type* or *scale/cross/type* *exclude/type* values = *types* *types* = list of types - *scale/cross/type* values = *shiftscale* *cmin* *wmin* + *scale/cross/type* values = *shiftscale* *cmin* *rmin* *shiftscale* = fraction of shifting in normal direction to preserve (unitless) *cmin* = minimum color function value required for scaling (unitless) - *wmin* = minimum local same-type support required for any shifting (unitless) + *rmin* = minimum local same-type weighted distance required for any shifting (unitless) *rho/sum* density evolution performed by a kernel summation values = none optional args = *self/mass* @@ -123,14 +123,14 @@ zero implies there is no shifting in the normal direction and a value of *scaleshift* of one implies no change in behavior. This scaling is only applied to atoms with a color function value greater than *cmin*. To handle scenarios of a small inclusion of one fluid type (e.g. a single atom) inside another, -the degree of same-type support is calculated +a weighted distance of same-type support is calculated .. math:: - W_{i,\mathrm{same}} = \sum_{j} W_{ij} \delta_{ij} + R_{i,\mathrm{same}} = \sum_{j} r_{ij} W_{ij} \delta_{ij} where :math:`\delta_{ij}` is zero if atoms :math:`i` and :math:`j` have different -types but unity otherwise. If :math:`W_{i,\mathrm{same}}` is ever less than the -specified value of *wmin*, shifting is turned off for particle :math:`i` +types but unity otherwise. If :math:`R_{i,\mathrm{same}}` is ever less than the +specified value of *rmin*, shifting is turned off for particle :math:`i` In systems with free surfaces (atom-vacuum), the *surface/detection* keyword can classify the location of particles as being within the bulk fluid, on a From 2dd6d30f072d804663d388e560812e0dbf8fd46d Mon Sep 17 00:00:00 2001 From: jtclemm Date: Fri, 14 Nov 2025 14:10:06 -0700 Subject: [PATCH 15/43] EdLimiting cv to be > 0 --- src/RHEO/fix_rheo_thermal.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index 61f4ade90ab..64b93b7f73e 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -128,7 +128,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error); double cv_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (cv_one < 0.0) error->all(FLERR, "The specific heat must be positive"); + if (cv_one <= 0.0) error->all(FLERR, "The specific heat must be greater than zero"); iarg += 2; for (i = nlo; i <= nhi; i++) { From 1582e48d552ec641f3e47ed6aa929a9d479e7222 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Sat, 15 Nov 2025 13:10:55 -0700 Subject: [PATCH 16/43] Fixed index bug in args, clarifying incrementing --- src/RHEO/fix_rheo.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index a9640bafe94..450b877ef75 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -126,16 +126,15 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : shift_scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp); shift_rmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp); - iarg += 3; + iarg += 4; } else if (strcmp(arg[iarg + 1], "exclude/type") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error); utils::bounds(FLERR, arg[iarg + 2], 1, n, nlo, nhi, error); for (i = nlo; i <= nhi; i++) shift_type[i] = 0; - iarg += 1; + iarg += 2; } else { break; } - iarg += 1; } } else if (strcmp(arg[iarg], "thermal") == 0) { thermal_flag = 1; @@ -146,25 +145,26 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : surface_style = COORDINATION; zmin_surface = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); zmin_splash = utils::inumeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 3; } else if (strcmp(arg[iarg + 1], "divergence") == 0) { surface_style = DIVR; divr_surface = utils::numeric(FLERR, arg[iarg + 2], false, lmp); zmin_splash = utils::inumeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 3; } else { error->all(FLERR, "Illegal surface/detection option in fix rheo, {}", arg[iarg + 1]); } - iarg += 3; } else if (strcmp(arg[iarg], "interface/reconstruct") == 0) { interface_flag = 1; } else if (strcmp(arg[iarg], "rho/sum") == 0) { rhosum_flag = 1; - while (iarg < narg) { // optional sub-arguments + while (iarg + 1 < narg) { // optional sub-arguments if (strcmp(arg[iarg + 1], "self/mass") == 0) { rhosum_self_mass_flag = 1; + iarg += 1; } else { break; } - iarg += 1; } } else if (strcmp(arg[iarg], "density") == 0) { if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo density", error); From 8a036a97f6c109c26eb120efbbef973ebaf110a3 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 17 Nov 2025 12:33:45 -0700 Subject: [PATCH 17/43] Avoiding int overflow in compute temp/sphere in large systems --- src/compute_temp_sphere.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/compute_temp_sphere.cpp b/src/compute_temp_sphere.cpp index 65c2cdf3cc1..33df1721506 100644 --- a/src/compute_temp_sphere.cpp +++ b/src/compute_temp_sphere.cpp @@ -126,7 +126,7 @@ void ComputeTempSphere::setup() void ComputeTempSphere::dof_compute() { - int count, count_all; + int count; adjust_dof_fix(); natoms_temp = group->count(igroup); @@ -168,7 +168,9 @@ void ComputeTempSphere::dof_compute() } } - MPI_Allreduce(&count, &count_all, 1, MPI_INT, MPI_SUM, world); + bigint count_single = count; + bigint count_all; + MPI_Allreduce(&count_single, &count_all, 1, MPI_LMP_BIGINT, MPI_SUM, world); dof = count_all; // additional adjustments to dof @@ -210,7 +212,8 @@ void ComputeTempSphere::dof_compute() } } - MPI_Allreduce(&count, &count_all, 1, MPI_INT, MPI_SUM, world); + count_single = count; + MPI_Allreduce(&count_single, &count_all, 1, MPI_LMP_BIGINT, MPI_SUM, world); dof -= count_all; } From ffdb39458d4a9d7ae4991c2d104e43c07a973ca9 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Thu, 18 Dec 2025 11:13:24 -0700 Subject: [PATCH 18/43] Adding access to the strain tensor --- src/EXTRA-FIX/fix_nonaffine_displacement.cpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp index 3329e0f05f6..243a6fec975 100644 --- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -146,6 +146,7 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char * comm_reverse = 0; comm_forward = 0; if (nad_style == D2MIN) { + size_peratom_cols = 9; comm_reverse = 18; comm_forward = 9; } @@ -195,7 +196,7 @@ void FixNonaffineDisplacement::post_constructor() grow_arrays(atom->nmax); for (int i = 0; i < atom->nlocal; i++) - for (int j = 0; j < 3; j++) array_atom[i][j] = 0.0; + for (int j = 0; j < size_peratom_cols; j++) array_atom[i][j] = 0.0; } /* ---------------------------------------------------------------------- */ @@ -627,6 +628,12 @@ void FixNonaffineDisplacement::calculate_D2Min() array_atom[i][0] = sqrt(D2min[i]); array_atom[i][1] = evol; array_atom[i][2] = edev; + array_atom[i][3] = E[0][0]; + array_atom[i][4] = E[1][1]; + array_atom[i][5] = E[2][2]; + array_atom[i][6] = E[0][1]; + array_atom[i][7] = E[0][2]; + array_atom[i][8] = E[1][2]; } } From 2136c1f97fba83bb07d1bb7d8752d90400e4aab6 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 6 Jan 2026 14:45:44 -0700 Subject: [PATCH 19/43] documenting d2min strain tensor output --- doc/src/fix_nonaffine_displacement.rst | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst index aaf8987e660..c46ec6ba839 100644 --- a/doc/src/fix_nonaffine_displacement.rst +++ b/doc/src/fix_nonaffine_displacement.rst @@ -112,14 +112,16 @@ The reference state is saved to :doc:`binary restart files `. None of the :doc:`fix_modify ` options are relevant to this fix. -This fix computes a peratom array with 3 columns, which can be accessed -by indices 1-3 using any command that uses per-atom values from a fix -as input. +This fix computes a peratom array with either 3 or 9 columns, which can +be accessed by indices 1-9 using any command that uses per-atom values +from a fix as input. For the *integrated* style, the three columns are the nonaffine displacements in the x, y, and z directions. For the *d2min* style, -the three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`, the -volumetric strain, and the deviatoric strain. +the first three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`, +the volumetric strain, and the deviatoric strain. The following 6 +columns are the xx, yy, zz, xy, xz, and yz components of the calculated +strain tensor. Restrictions """""""""""" From d4d3166508a9b01cfa4de14ee01a2cacd00f8fa7 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 6 Jan 2026 14:45:57 -0700 Subject: [PATCH 20/43] Adding back in dropped rheo references --- doc/src/set.rst | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/doc/src/set.rst b/doc/src/set.rst index 10a6c0993d7..5a6995ee09b 100644 --- a/doc/src/set.rst +++ b/doc/src/set.rst @@ -27,12 +27,12 @@ Syntax or *density* or *density/disc* or *diameter* or *dihedral* or *dipole* or *dipole/random* or *dpd/theta* or *edpd/cv* or *edpd/temp* or *epsilon* or *image* or *improper* or *length* or *mass* or *mol* or - *omega* or *quat* or *quat/random* or *radius/electron* or *shape* or - *smd/contact/radius* or *smd/mass/density* or *sph/cv* or *sph/e* or - *sph/rho* or *spin/atom* or *spin/atom/random* or *spin/electron* or - *temperature* or *theta* or *theta/random* or *tri* or *type* or - *type/fraction* or *type/ratio* or *type/subset* or *volume* or *vx* - or *vy* or *vz* or *x* or *y* or *z* or *i_name* or *d_name* or + *omega* or *quat* or *quat/random* or *radius/electron* or *rheo/rho* or + *rheo/status* or *shape* or *smd/contact/radius* or *smd/mass/density* or + *sph/cv* or *sph/e* or *sph/rho* or *spin/atom* or *spin/atom/random* or + *spin/electron* or *temperature* or *theta* or *theta/random* or *tri* or + *type* or *type/fraction* or *type/ratio* or *type/subset* or *volume* or + *vx* or *vy* or *vz* or *x* or *y* or *z* or *i_name* or *d_name* or *i2_name* or *d2_name* .. parsed-literal:: @@ -97,6 +97,10 @@ Syntax *radius/electron* value = eradius eradius = electron radius (or fixed-core radius) (distance units) value can be an atom-style variable (see below) + *rheo/rho* value = density of RHEO particles (mass/distance\^3) + value can be an atom-style variable (see below) + *rheo/status* value = status or phase of RHEO particles (unitless) + value can be an atom-style variable (see below) *shape* values = Sx Sy Sz Sx,Sy,Sz = 3 diameters of ellipsoid (distance units) any of Sx,Sy,Sz can be an atom-style variable (see below) @@ -512,6 +516,10 @@ use of an atom-style variable. Keyword *radius/electron* uses the specified value to set the radius of electrons or fixed cores. +Keywords *rheo/rho* and *rheo/status* set the density and the status of +rheo particles. In particular, one can only set the phase in the status +as described by the :doc:`RHEO howto page `. + Keyword *shape* sets the size and shape of the selected atoms. The particles must be ellipsoids as defined by the :doc:`atom_style ellipsoid ` command. The *Sx*, *Sy*, *Sz* settings are From 55b23346cbf1b9e82cef60f881da1b7a13ea65d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 11:47:30 +0100 Subject: [PATCH 21/43] Add files via upload Example of using the Streitz-Mintmire model for GaN with smooth termination of charge transfer at the cutoff. --- examples/streitz/GaN.sm | 3 ++ examples/streitz/GaN.sm+tersoff.mod | 28 +++++++++++ examples/streitz/in.gan | 73 +++++++++++++++++++++++++++++ examples/streitz/scmin.mod | 37 +++++++++++++++ 4 files changed, 141 insertions(+) create mode 100644 examples/streitz/GaN.sm create mode 100644 examples/streitz/GaN.sm+tersoff.mod create mode 100644 examples/streitz/in.gan create mode 100644 examples/streitz/scmin.mod diff --git a/examples/streitz/GaN.sm b/examples/streitz/GaN.sm new file mode 100644 index 00000000000..d5eb6cc929e --- /dev/null +++ b/examples/streitz/GaN.sm @@ -0,0 +1,3 @@ +# chi (eV), J (eV), gamma (1/ang), zeta (1/ang), Z (e) +Ga -1.25 14.51 0 1.2333912356832242 0 +N 3.38 6.91 0 0.8509319105529678 0 diff --git a/examples/streitz/GaN.sm+tersoff.mod b/examples/streitz/GaN.sm+tersoff.mod new file mode 100644 index 00000000000..0fbc7336116 --- /dev/null +++ b/examples/streitz/GaN.sm+tersoff.mod @@ -0,0 +1,28 @@ +# DATE: 2026-01-27 CONTRIBUTOR: Roman Groger +# CITATION: R. Groger, J. Fikar, Acta Mater. XX, XXXX (2026) +# +# Modified Tersoff potential parameters for GaN describing short-range interactions in ionic-covalent +# model of GaN developed within the framework of the Streitz-Mintmire formalism. This file must be +# combined with the electrostatic energy to form a hybrid potential. The electrostatic parameters +# of Ga and N are stored in the file GaN.sm. +# +# The parameters of homobonds are taken from Nord et al., J. Phys. Cond. Mat 15, 5649 (2003). +# Fitting of Ga-N parameters and the development of the ionic-covalent model for GaN is described in +# Groger & Fikar, Acta Mater. XX, XXXX (2026) +# +# Multiple entries can be added to this file, LAMMPS reads the ones it needs. +# These entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless +# +# Format of a single entry (one or more lines): +# element 1, element 2, element 3, beta, alpha, h, eta, beta_ters, lambda2, B, R, D, lambda1, A, n, c1, c2, c3, c4, c5 + +Ga Ga Ga 1.0 1.846 -0.3013 1.0 1.0 1.4497 410.132 2.87 0.15 1.60916 535.199 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 +N N N 1.0 0.0 -0.045238 1.0 1.0 2.38426 423.769 2.2 0.2 3.55779 1044.77 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 +Ga Ga N 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +Ga N N 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +N Ga Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +N Ga N 1.0 0.0 -0.045238 1.0 1.0 0.0 0.0 2.2 0.2 0.0 0.0 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 +N N Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +Ga N Ga 1.0 1.846 -0.3013 1.0 1.0 0.0 0.0 2.87 0.15 0.0 0.0 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 diff --git a/examples/streitz/in.gan b/examples/streitz/in.gan new file mode 100644 index 00000000000..111bb6a0237 --- /dev/null +++ b/examples/streitz/in.gan @@ -0,0 +1,73 @@ +variable a0 equal 3.0 # trial values +variable c0 equal 5.0 +variable u0 equal 0.38 +variable Efu_inf equal -0.9183986628325442 # energy of a pair of isolated Ga and N ions + +units metal +dimension 3 +boundary p p p +atom_style charge +atom_modify map yes + +lattice custom ${a0} & + a1 1.0 0.0 0.0 & + a2 -0.5 0.866025403784439 0.0 & + a3 0.0 0.0 $(v_c0/v_a0) & + spacing 1 0.866025403784439 $(v_c0/v_a0) & + basis 0.33333333 0.66666667 0 & + basis 0.66666667 0.33333333 0.5 & + basis 0.33333333 0.66666667 ${u0} & + basis 0.66666667 0.33333333 $(v_u0+0.5) +region box prism 0 1 0 1 0 1 -0.5 0 0 units lattice +create_box 2 box +create_atoms 2 box basis 1 1 basis 2 1 basis 3 2 basis 4 2 +mass * 1.0 + +set type 1 charge 1.0 +set type 2 charge -1.0 +group g1 type 1 +group g2 type 2 + +compute eat all pe/atom +compute q1all g1 property/atom q +compute q1 g1 reduce ave c_q1all + +pair_style hybrid/overlay tersoff/mod coul/streitz 12.0 dsf 0.3 2.0 +pair_coeff * * tersoff/mod GaN.sm+tersoff.mod Ga N +pair_coeff * * coul/streitz GaN.sm Ga N + +fix fq all qeq/slater 1 12.0 1e-6 1000 coul/streitz alpha 0.3 dsf 2.0 + +fix 1 all box/relax aniso 0.0 + +#dump 1 all cfg 100 dump_*.cfg mass type xs ys zs c_eat q +#dump_modify 1 element Ga N + +thermo 100 +thermo_style custom step pe fmax press vol lx ly lz +include scmin.mod # self-consistent minimization & charge equilibration + +variable Efu_nocorr equal 2*pe/count(all) +variable Efu equal ${Efu_nocorr}-${Efu_inf} +variable a equal (lx+ly*2/sqrt(3))/2 +variable c equal lz +variable p equal press/1e4 +variable q1 equal c_q1 + +print '--------------------------------------------------------------------------------------------' +print 'Energy (no correction): E = ${Efu_nocorr} eV/fu' +print 'Energy: E-Einf = ${Efu} eV/fu' +print 'Pressure: p = $p GPa' +print 'Energy of a pair of isolated ions: Einf = ${Efu_inf} eV/fu' +print 'Lattice parameters: a = $a A' +print ' c = $c A' +print 'Cationic charge: q = ${q1}*e' +print '--------------------------------------------------------------------------------------------' + +#include polariz.mod +#include polariz_diff.mod +#include polariz_wz.mod # no relaxation, fixed cell shape +#include polariz_wz_relax.mod # with relaxation - transforms into different structures? + + + diff --git a/examples/streitz/scmin.mod b/examples/streitz/scmin.mod new file mode 100644 index 00000000000..00279c5f468 --- /dev/null +++ b/examples/streitz/scmin.mod @@ -0,0 +1,37 @@ +# QEq cannot be run together with minimize, becase each new calculation of charges alters the energy +# hypersurface on which we are looking for a minimum. Minimization frequently gets stuck with the +# message "linesearch alpha is zero" before reaching zero pressure. +# +# Below is a self-consistency loop that minimizes the potential energy while distributing atomic +# charges, which replaces the minimize command. It works as follows: +# 1. For given atomic positions, run QEq to compute charges. The initial charges are q0 and the +# final ones q. +# 2. Fix charges (q) and run energy minimization with changing box shape to get new atomic +# coordinates. This gives zero pressure. +# 3. If max|q-q0| >= 0.01, set q0 := q and go to 1. +# If max|q-q0| < 0.01, we have equilibrium atomic positions and charges => exit. +# +# Local variables end with _ to avoid interference with the main script. + +variable dq_ atom abs(q-f_q0_) +compute dqmax_ all reduce max v_dq_ +variable dqmax_ equal c_dqmax_ + +variable iter loop 100 # number of iterations to equilibrate charges & get zero pressure +label self_consistency + fix q0_ all store/state 0 q + fix fq all qeq/slater 1 12.0 1e-6 1000 coul/streitz alpha 0.3 dsf 2.0 + run 0 # calculate charges -> nonzero pressure + + unfix fq + minimize 0 1e-6 10000 10000 # fix charges -> zero pressure + + print "***** self-consistent loop: iter=${iter}, dqmax=${dqmax_} *****" + if "${dqmax_} < 0.01" then "jump SELF break" + next iter +jump SELF self_consistency + +label break +unfix q0_ +uncompute dqmax_ +variable iter delete From 674fe7824335409e1a4f0391e19d07c4c4357d38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 11:54:49 +0100 Subject: [PATCH 22/43] Add files via upload tersoff/mod potential for GaN with charge equilibration using the Streitz-Mintmire model --- potentials/GaN.streitz | 3 +++ potentials/GaN.streitz+tersoff.mod | 28 ++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) create mode 100644 potentials/GaN.streitz create mode 100644 potentials/GaN.streitz+tersoff.mod diff --git a/potentials/GaN.streitz b/potentials/GaN.streitz new file mode 100644 index 00000000000..d5eb6cc929e --- /dev/null +++ b/potentials/GaN.streitz @@ -0,0 +1,3 @@ +# chi (eV), J (eV), gamma (1/ang), zeta (1/ang), Z (e) +Ga -1.25 14.51 0 1.2333912356832242 0 +N 3.38 6.91 0 0.8509319105529678 0 diff --git a/potentials/GaN.streitz+tersoff.mod b/potentials/GaN.streitz+tersoff.mod new file mode 100644 index 00000000000..304ba73ceee --- /dev/null +++ b/potentials/GaN.streitz+tersoff.mod @@ -0,0 +1,28 @@ +# DATE: 2026-01-27 CONTRIBUTOR: Roman Groger +# CITATION: R. Groger, J. Fikar, Acta Mater. XX, XXXX (2026) +# +# Modified Tersoff potential parameters for GaN describing short-range interactions in ionic-covalent +# model of GaN developed within the framework of the Streitz-Mintmire formalism. This file must be +# combined with the electrostatic energy to form a hybrid potential. The electrostatic parameters +# of Ga and N are stored in the file GaN.streitz. +# +# The parameters of homobonds are taken from Nord et al., J. Phys. Cond. Mat 15, 5649 (2003). +# Fitting of Ga-N parameters and the development of the ionic-covalent model for GaN is described in +# Groger & Fikar, Acta Mater. XX, XXXX (2026) +# +# Multiple entries can be added to this file, LAMMPS reads the ones it needs. +# These entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless +# +# Format of a single entry (one or more lines): +# element 1, element 2, element 3, beta, alpha, h, eta, beta_ters, lambda2, B, R, D, lambda1, A, n, c1, c2, c3, c4, c5 + +Ga Ga Ga 1.0 1.846 -0.3013 1.0 1.0 1.4497 410.132 2.87 0.15 1.60916 535.199 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 +N N N 1.0 0.0 -0.045238 1.0 1.0 2.38426 423.769 2.2 0.2 3.55779 1044.77 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 +Ga Ga N 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +Ga N N 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +N Ga Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +N Ga N 1.0 0.0 -0.045238 1.0 1.0 0.0 0.0 2.2 0.2 0.0 0.0 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 +N N Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +Ga N Ga 1.0 1.846 -0.3013 1.0 1.0 0.0 0.0 2.87 0.15 0.0 0.0 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 From 8ee1c2c9611778fecf22d009439a2093a1877c51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 11:56:04 +0100 Subject: [PATCH 23/43] Add files via upload --- examples/streitz/GaN.streitz | 3 +++ examples/streitz/GaN.streitz+tersoff.mod | 28 ++++++++++++++++++++++++ examples/streitz/in.gan | 4 ++-- 3 files changed, 33 insertions(+), 2 deletions(-) create mode 100644 examples/streitz/GaN.streitz create mode 100644 examples/streitz/GaN.streitz+tersoff.mod diff --git a/examples/streitz/GaN.streitz b/examples/streitz/GaN.streitz new file mode 100644 index 00000000000..d5eb6cc929e --- /dev/null +++ b/examples/streitz/GaN.streitz @@ -0,0 +1,3 @@ +# chi (eV), J (eV), gamma (1/ang), zeta (1/ang), Z (e) +Ga -1.25 14.51 0 1.2333912356832242 0 +N 3.38 6.91 0 0.8509319105529678 0 diff --git a/examples/streitz/GaN.streitz+tersoff.mod b/examples/streitz/GaN.streitz+tersoff.mod new file mode 100644 index 00000000000..304ba73ceee --- /dev/null +++ b/examples/streitz/GaN.streitz+tersoff.mod @@ -0,0 +1,28 @@ +# DATE: 2026-01-27 CONTRIBUTOR: Roman Groger +# CITATION: R. Groger, J. Fikar, Acta Mater. XX, XXXX (2026) +# +# Modified Tersoff potential parameters for GaN describing short-range interactions in ionic-covalent +# model of GaN developed within the framework of the Streitz-Mintmire formalism. This file must be +# combined with the electrostatic energy to form a hybrid potential. The electrostatic parameters +# of Ga and N are stored in the file GaN.streitz. +# +# The parameters of homobonds are taken from Nord et al., J. Phys. Cond. Mat 15, 5649 (2003). +# Fitting of Ga-N parameters and the development of the ionic-covalent model for GaN is described in +# Groger & Fikar, Acta Mater. XX, XXXX (2026) +# +# Multiple entries can be added to this file, LAMMPS reads the ones it needs. +# These entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless +# +# Format of a single entry (one or more lines): +# element 1, element 2, element 3, beta, alpha, h, eta, beta_ters, lambda2, B, R, D, lambda1, A, n, c1, c2, c3, c4, c5 + +Ga Ga Ga 1.0 1.846 -0.3013 1.0 1.0 1.4497 410.132 2.87 0.15 1.60916 535.199 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 +N N N 1.0 0.0 -0.045238 1.0 1.0 2.38426 423.769 2.2 0.2 3.55779 1044.77 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 +Ga Ga N 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +Ga N N 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +N Ga Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +N Ga N 1.0 0.0 -0.045238 1.0 1.0 0.0 0.0 2.2 0.2 0.0 0.0 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 +N N Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 +Ga N Ga 1.0 1.846 -0.3013 1.0 1.0 0.0 0.0 2.87 0.15 0.0 0.0 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 diff --git a/examples/streitz/in.gan b/examples/streitz/in.gan index 111bb6a0237..8e6a477658b 100644 --- a/examples/streitz/in.gan +++ b/examples/streitz/in.gan @@ -33,8 +33,8 @@ compute q1all g1 property/atom q compute q1 g1 reduce ave c_q1all pair_style hybrid/overlay tersoff/mod coul/streitz 12.0 dsf 0.3 2.0 -pair_coeff * * tersoff/mod GaN.sm+tersoff.mod Ga N -pair_coeff * * coul/streitz GaN.sm Ga N +pair_coeff * * tersoff/mod GaN.streitz+tersoff.mod Ga N +pair_coeff * * coul/streitz GaN.streitz Ga N fix fq all qeq/slater 1 12.0 1e-6 1000 coul/streitz alpha 0.3 dsf 2.0 From 1fc188ec5a473928911987c32038f06bca6082d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 12:00:37 +0100 Subject: [PATCH 24/43] Delete examples/streitz/GaN.sm --- examples/streitz/GaN.sm | 3 --- 1 file changed, 3 deletions(-) delete mode 100644 examples/streitz/GaN.sm diff --git a/examples/streitz/GaN.sm b/examples/streitz/GaN.sm deleted file mode 100644 index d5eb6cc929e..00000000000 --- a/examples/streitz/GaN.sm +++ /dev/null @@ -1,3 +0,0 @@ -# chi (eV), J (eV), gamma (1/ang), zeta (1/ang), Z (e) -Ga -1.25 14.51 0 1.2333912356832242 0 -N 3.38 6.91 0 0.8509319105529678 0 From 06588a30f499f6cfd53bb21927138bc30cc7eb02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 12:00:48 +0100 Subject: [PATCH 25/43] Delete examples/streitz/GaN.sm+tersoff.mod --- examples/streitz/GaN.sm+tersoff.mod | 28 ---------------------------- 1 file changed, 28 deletions(-) delete mode 100644 examples/streitz/GaN.sm+tersoff.mod diff --git a/examples/streitz/GaN.sm+tersoff.mod b/examples/streitz/GaN.sm+tersoff.mod deleted file mode 100644 index 0fbc7336116..00000000000 --- a/examples/streitz/GaN.sm+tersoff.mod +++ /dev/null @@ -1,28 +0,0 @@ -# DATE: 2026-01-27 CONTRIBUTOR: Roman Groger -# CITATION: R. Groger, J. Fikar, Acta Mater. XX, XXXX (2026) -# -# Modified Tersoff potential parameters for GaN describing short-range interactions in ionic-covalent -# model of GaN developed within the framework of the Streitz-Mintmire formalism. This file must be -# combined with the electrostatic energy to form a hybrid potential. The electrostatic parameters -# of Ga and N are stored in the file GaN.sm. -# -# The parameters of homobonds are taken from Nord et al., J. Phys. Cond. Mat 15, 5649 (2003). -# Fitting of Ga-N parameters and the development of the ionic-covalent model for GaN is described in -# Groger & Fikar, Acta Mater. XX, XXXX (2026) -# -# Multiple entries can be added to this file, LAMMPS reads the ones it needs. -# These entries are in LAMMPS "metal" units: -# A,B = eV; lambda1,lambda2 = 1/Angstroms; R,D = Angstroms -# other quantities are unitless -# -# Format of a single entry (one or more lines): -# element 1, element 2, element 3, beta, alpha, h, eta, beta_ters, lambda2, B, R, D, lambda1, A, n, c1, c2, c3, c4, c5 - -Ga Ga Ga 1.0 1.846 -0.3013 1.0 1.0 1.4497 410.132 2.87 0.15 1.60916 535.199 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 -N N N 1.0 0.0 -0.045238 1.0 1.0 2.38426 423.769 2.2 0.2 3.55779 1044.77 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 -Ga Ga N 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -Ga N N 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -N Ga Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -N Ga N 1.0 0.0 -0.045238 1.0 1.0 0.0 0.0 2.2 0.2 0.0 0.0 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 -N N Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -Ga N Ga 1.0 1.846 -0.3013 1.0 1.0 0.0 0.0 2.87 0.15 0.0 0.0 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 From b51bb5f5cb1ca1b75afbe151de7a5f16a0944d52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 13:41:40 +0100 Subject: [PATCH 26/43] Update fix_qeq.rst Documentation on new options for specifying taper widths to smooth terminate Coulomb integrals when calculating electrostatic interactions. --- doc/src/fix_qeq.rst | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_qeq.rst b/doc/src/fix_qeq.rst index ac0bca84f16..42bc9d3c946 100644 --- a/doc/src/fix_qeq.rst +++ b/doc/src/fix_qeq.rst @@ -42,7 +42,9 @@ Syntax .. parsed-literal:: - *alpha* value = Slater type orbital exponent (qeq/slater only) + *alpha* value = Slater type orbital exponent (qeq/slater only). Can be followed by optional arguments: + *wolf* value = width of taper to terminate Coulomb integrals for the Wolf summation (default value is zero) + *dsf* value = width of taper to terminate Coulomb integrals for the Fennell-Gezelter summation (default value is zero) *cdamp* value = damping parameter for Coulomb interactions (qeq/ctip only) *maxrepeat* value = number of equilibration cycles allowed to ensure no atoms cross charge bounds (qeq/ctip only) *qdamp* value = damping factor for damped dynamics charge solver (qeq/dynamic and qeq/fire only) @@ -57,6 +59,10 @@ Examples fix 1 all qeq/point 1 10 1.0e-6 200 param.qeq1 fix 1 qeq qeq/shielded 1 8 1.0e-6 100 param.qeq2 fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 + fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 wolf + fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 wolf 2.0 + fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 dsf + fix 1 all qeq/slater 5 10 1.0e-6 100 params alpha 0.2 dsf 2.0 fix 1 all qeq/ctip 1 12 1.0e-8 100 coul/ctip cdamp 0.30 maxrepeat 10 fix 1 qeq qeq/dynamic 1 12 1.0e-3 100 my_qeq fix 1 all qeq/fire 1 10 1.0e-3 100 my_qeq qdamp 0.2 qstep 0.1 @@ -249,6 +255,8 @@ larger sizes, and *qeq/fire* is faster than *qeq/dynamic*\ . arbitrary choices of these parameters. We do not develop these QEq parameters. See the examples/qeq directory for some examples. +In the previous versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing 1/r using a damped potential erfc(alpha*r)/r with the parameter *alpha* controlling the rate of decay. However, any choice of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potentials due to Wolf et al. (*wolf*) and its extension by Fennell and Gezelter (*dsf*) solve this problem, but they were not implemented for charge equilibration in previous versions of LAMMPS. An extension was implemented to specify the width of taper (see Murty et al.) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of the taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. + Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From 4a2d8f9b49036cbe90b29cab58568f3645739c57 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 14:01:56 +0100 Subject: [PATCH 27/43] Update fix_qeq.rst Updated documentation - usage of optional arguments wolf and dsf and user-defined taper width to smoothly terminate Coulomb integrals. --- doc/src/fix_qeq.rst | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_qeq.rst b/doc/src/fix_qeq.rst index 42bc9d3c946..2b86b3017e6 100644 --- a/doc/src/fix_qeq.rst +++ b/doc/src/fix_qeq.rst @@ -43,8 +43,8 @@ Syntax .. parsed-literal:: *alpha* value = Slater type orbital exponent (qeq/slater only). Can be followed by optional arguments: - *wolf* value = width of taper to terminate Coulomb integrals for the Wolf summation (default value is zero) - *dsf* value = width of taper to terminate Coulomb integrals for the Fennell-Gezelter summation (default value is zero) + *wolf* value = width of taper to terminate Coulomb integrals for the Wolf summation (default value is zero) + *dsf* value = width of taper to terminate Coulomb integrals for the Fennell-Gezelter summation (default value is zero) *cdamp* value = damping parameter for Coulomb interactions (qeq/ctip only) *maxrepeat* value = number of equilibration cycles allowed to ensure no atoms cross charge bounds (qeq/ctip only) *qdamp* value = damping factor for damped dynamics charge solver (qeq/dynamic and qeq/fire only) @@ -255,7 +255,7 @@ larger sizes, and *qeq/fire* is faster than *qeq/dynamic*\ . arbitrary choices of these parameters. We do not develop these QEq parameters. See the examples/qeq directory for some examples. -In the previous versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing 1/r using a damped potential erfc(alpha*r)/r with the parameter *alpha* controlling the rate of decay. However, any choice of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potentials due to Wolf et al. (*wolf*) and its extension by Fennell and Gezelter (*dsf*) solve this problem, but they were not implemented for charge equilibration in previous versions of LAMMPS. An extension was implemented to specify the width of taper (see Murty et al.) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of the taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. +In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -322,3 +322,15 @@ Physical Chemistry, 105, 9396-9049 (2001) .. _Shan: **(QEq/Fire)** T.-R. Shan, A. P. Thompson, S. J. Plimpton, in preparation + +.. _Wolf: + +**(Wolf)** D. Wolf, P. Keblinski, S. R. Phillpot, J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999). + +.. _Fennell: + +**(Fennell)** J. Fennell, J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006). + +.. _Mei: + +**(Mei)** J. Mei, J. W. Davenport, G. W. Fernando, Phys. Rev. B 43, 4653 (1991). From 8a9cbaa118ab51d85e775b01b566c022f7e17df7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 14:08:38 +0100 Subject: [PATCH 28/43] Update pair_coul.rst --- doc/src/pair_coul.rst | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index c6735612379..2eeb1a3dfc8 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -146,6 +146,11 @@ Examples pair_style coul/streitz 12.0 wolf 0.30 pair_coeff * * AlO.streitz Al O + pair_style coul/streitz 12.0 wolf 0.3 2.0 + pair_style coul/streitz 12.0 dsf 0.3 + pair_style coul/streitz 12.0 dsf 0.3 2.0 + pair_coeff * * GaN.streitz Ga N + pair_style tip4p/cut 1 2 7 8 0.15 12.0 pair_coeff * * @@ -302,7 +307,7 @@ command doc page. Alternatively *qfile* can be replaced by "coul/streitz", in which case the fix will extract QEq parameters from the coul/streitz pair style itself. -See the examples/strietz directory for an example input script that +See the examples/streitz directory for an example input script that uses the Streitz-Mintmire potential. The potentials directory has the AlO.eam.alloy and AlO.streitz potential files used by the example. From 07b3ffde22207d9d522ea7b16308aa91df03baec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 14:41:10 +0100 Subject: [PATCH 29/43] Update pair_coul.rst Extension of pair coul/streitz to include smooth taper at the cutoff, GaN as an example. --- doc/src/pair_coul.rst | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index 2eeb1a3dfc8..1d27e4c5309 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -146,8 +146,9 @@ Examples pair_style coul/streitz 12.0 wolf 0.30 pair_coeff * * AlO.streitz Al O + pair_style coul/streitz 12.0 wolf 0.3 # default taper width is zero pair_style coul/streitz 12.0 wolf 0.3 2.0 - pair_style coul/streitz 12.0 dsf 0.3 + pair_style coul/streitz 12.0 dsf 0.3 # default taper width is zero pair_style coul/streitz 12.0 dsf 0.3 2.0 pair_coeff * * GaN.streitz Ga N @@ -311,13 +312,25 @@ See the examples/streitz directory for an example input script that uses the Streitz-Mintmire potential. The potentials directory has the AlO.eam.alloy and AlO.streitz potential files used by the example. +In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. + Note that the Streiz-Mintmire potential is generally used for oxides, but there is no conceptual problem with extending it to nitrides and -carbides (such as SiC, TiN). Pair coul/strietz used by itself or with +carbides (such as SiC, TiN). Pair coul/streitz used by itself or with any other pair style such as EAM, MEAM, Tersoff, or LJ in hybrid/overlay mode. To do this, you would need to provide a Streitz-Mintmire parameterization for the material being modeled. +An implementation of the Streitz-Mintmire potential for GaN due to :ref:`(Groger and Fikar) ` can be found in the examples/streitz directory. The electrostatic parameters of Ga and N are stored in file GaN.streitz and the short-range tersoff/mod potential in the file GaN.streitz+tersoff.mod. The total potential must be specified as: + +.. code-block:: LAMMPS + + pair_style hybrid/overlay tersoff/mod coul/streitz 12.0 dsf 0.3 2.0 + pair_coeff * * tersoff/mod GaN.streitz+tersoff.mod Ga N + pair_coeff * * coul/streitz GaN.streitz Ga N + +where the last three parameters specify the method of calculation of Coulomb interactions (*ewald*, *wolf* or *dsf*), the value of *alpha*, and the width of taper to smoothly terminate the Coulomb integrals. + ---------- Pair style *coul/exclude* computes Coulombic interactions like *coul/cut* @@ -482,3 +495,11 @@ J Chemical Physics, 162, 054709 (2025). **(Jorgensen)** Jorgensen, Chandrasekhar, Madura, Impey, Klein, J Chem Phys, 79, 926 (1983). + +.. _Mei1: + +**(Mei)** J. Mei, J. W. Davenport, G. W. Fernando, Phys. Rev. B 43, 4653 (1991). + +.. _Groger1: + +**(Groger)** R. Groger, J. Fikar, Acta Mater. (2026) in review. From 1d632278995d392a70272a27e0b9f00ef0b86354 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 14:42:43 +0100 Subject: [PATCH 30/43] Update pair_coul.rst --- doc/src/pair_coul.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index 1d27e4c5309..511713279d9 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -321,7 +321,7 @@ any other pair style such as EAM, MEAM, Tersoff, or LJ in hybrid/overlay mode. To do this, you would need to provide a Streitz-Mintmire parameterization for the material being modeled. -An implementation of the Streitz-Mintmire potential for GaN due to :ref:`(Groger and Fikar) ` can be found in the examples/streitz directory. The electrostatic parameters of Ga and N are stored in file GaN.streitz and the short-range tersoff/mod potential in the file GaN.streitz+tersoff.mod. The total potential must be specified as: +An implementation of the Streitz-Mintmire potential for GaN due to :ref:`(Groger and Fikar) ` can be found in the examples/streitz directory. The electrostatic parameters of Ga and N are stored in file GaN.streitz and the short-range tersoff/mod potential in the file GaN.streitz+tersoff.mod. The total potential must be specified as: .. code-block:: LAMMPS From 9d62a7bf66fc66b196cda84e0f800d59eef9336b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 14:43:09 +0100 Subject: [PATCH 31/43] Update fix_qeq.rst --- doc/src/fix_qeq.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_qeq.rst b/doc/src/fix_qeq.rst index 2b86b3017e6..891098721c1 100644 --- a/doc/src/fix_qeq.rst +++ b/doc/src/fix_qeq.rst @@ -255,7 +255,7 @@ larger sizes, and *qeq/fire* is faster than *qeq/dynamic*\ . arbitrary choices of these parameters. We do not develop these QEq parameters. See the examples/qeq directory for some examples. -In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. +In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From efb5e07d73d373cc5fe88a801db2ea0a6cb63c80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 17:17:07 +0100 Subject: [PATCH 32/43] Update doc/src/pair_coul.rst Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- doc/src/pair_coul.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index 511713279d9..ca29b2e6645 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -314,7 +314,7 @@ AlO.eam.alloy and AlO.streitz potential files used by the example. In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. -Note that the Streiz-Mintmire potential is generally used for oxides, +Note that the Streitz-Mintmire potential is generally used for oxides, but there is no conceptual problem with extending it to nitrides and carbides (such as SiC, TiN). Pair coul/streitz used by itself or with any other pair style such as EAM, MEAM, Tersoff, or LJ in From 1e7354d56cbe10517b6ae77f02a5474696b27251 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 17:17:39 +0100 Subject: [PATCH 33/43] Update doc/src/pair_coul.rst Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- doc/src/pair_coul.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index ca29b2e6645..eb7bbf36e64 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -312,7 +312,7 @@ See the examples/streitz directory for an example input script that uses the Streitz-Mintmire potential. The potentials directory has the AlO.eam.alloy and AlO.streitz potential files used by the example. -In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. +In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. Note that the Streitz-Mintmire potential is generally used for oxides, but there is no conceptual problem with extending it to nitrides and From 050516bfffcbcf64004e6832f7f007373c3addb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:00:56 +0100 Subject: [PATCH 34/43] Update fix_qeq.rst --- doc/src/fix_qeq.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/src/fix_qeq.rst b/doc/src/fix_qeq.rst index 891098721c1..ebc42d43885 100644 --- a/doc/src/fix_qeq.rst +++ b/doc/src/fix_qeq.rst @@ -255,7 +255,9 @@ larger sizes, and *qeq/fire* is faster than *qeq/dynamic*\ . arbitrary choices of these parameters. We do not develop these QEq parameters. See the examples/qeq directory for some examples. -In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. +.. versionadded:: TBD + +In previous versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" From 2a0ae755d17ef785bcb41dc72f02faa4335b70c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:02:06 +0100 Subject: [PATCH 35/43] Update scmin.mod --- examples/streitz/scmin.mod | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/streitz/scmin.mod b/examples/streitz/scmin.mod index 00279c5f468..50e504aadd4 100644 --- a/examples/streitz/scmin.mod +++ b/examples/streitz/scmin.mod @@ -1,4 +1,4 @@ -# QEq cannot be run together with minimize, becase each new calculation of charges alters the energy +# QEq cannot be run together with minimize, because each new calculation of charges alters the energy # hypersurface on which we are looking for a minimum. Minimization frequently gets stuck with the # message "linesearch alpha is zero" before reaching zero pressure. # From ccf6e8a53d7dcf33232b32e7ba8e6d87b9854f85 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:05:17 +0100 Subject: [PATCH 36/43] Update pair_coul_streitz.cpp --- src/KSPACE/pair_coul_streitz.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/KSPACE/pair_coul_streitz.cpp b/src/KSPACE/pair_coul_streitz.cpp index 9206c41796c..ed958728741 100644 --- a/src/KSPACE/pair_coul_streitz.cpp +++ b/src/KSPACE/pair_coul_streitz.cpp @@ -50,6 +50,7 @@ PairCoulStreitz::PairCoulStreitz(LAMMPS *lmp) : Pair(lmp) one_coeff = 1; nmax = 0; drtap = 0.0; + dsfflag = 0; params = nullptr; } From 3a9f92d4464127e89a83a787497164239bbd144a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:09:05 +0100 Subject: [PATCH 37/43] Update pair_coul.rst --- doc/src/pair_coul.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index eb7bbf36e64..205af881571 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -312,8 +312,6 @@ See the examples/streitz directory for an example input script that uses the Streitz-Mintmire potential. The potentials directory has the AlO.eam.alloy and AlO.streitz potential files used by the example. -In older versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. - Note that the Streitz-Mintmire potential is generally used for oxides, but there is no conceptual problem with extending it to nitrides and carbides (such as SiC, TiN). Pair coul/streitz used by itself or with @@ -321,6 +319,10 @@ any other pair style such as EAM, MEAM, Tersoff, or LJ in hybrid/overlay mode. To do this, you would need to provide a Streitz-Mintmire parameterization for the material being modeled. +.. versionchanged:: TBD + +In previous versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. + An implementation of the Streitz-Mintmire potential for GaN due to :ref:`(Groger and Fikar) ` can be found in the examples/streitz directory. The electrostatic parameters of Ga and N are stored in file GaN.streitz and the short-range tersoff/mod potential in the file GaN.streitz+tersoff.mod. The total potential must be specified as: .. code-block:: LAMMPS From 283e142a5cf63a42dba971ae788d0f3f7d4759d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:12:52 +0100 Subject: [PATCH 38/43] Update src/KSPACE/pair_coul_streitz.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/KSPACE/pair_coul_streitz.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/KSPACE/pair_coul_streitz.cpp b/src/KSPACE/pair_coul_streitz.cpp index ed958728741..4b22920e16c 100644 --- a/src/KSPACE/pair_coul_streitz.cpp +++ b/src/KSPACE/pair_coul_streitz.cpp @@ -51,6 +51,7 @@ PairCoulStreitz::PairCoulStreitz(LAMMPS *lmp) : Pair(lmp) nmax = 0; drtap = 0.0; dsfflag = 0; + dsfflag = 0; params = nullptr; } From 41558eae6213a84abd61ad0ee1a72588545a9a6c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:13:32 +0100 Subject: [PATCH 39/43] Update src/QEQ/fix_qeq_slater.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/QEQ/fix_qeq_slater.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index be0e2346ab3..c63e5045540 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -76,7 +76,6 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg iarg += 2; } else { error->all(FLERR, "Unknown fix qeq/slater keyword: {}", arg[iarg]); - continue; } } From 058e9e495233fcd8122f8da54a76f46396c341b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:19:26 +0100 Subject: [PATCH 40/43] Update fix_qeq_slater.cpp --- src/QEQ/fix_qeq_slater.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index c63e5045540..c9dd0da19f7 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -57,14 +57,16 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg if (strcmp(arg[iarg],"wolf") == 0) { // type of potential (0=erfc(r)/r, 1=Wolf, 2=Fennell & Gezelter) vtype = 1; iarg++; - if (iarg < narg) { + // optional taper width; do not treat following keyword (e.g. "warn") as numeric + if (iarg < narg && strcmp(arg[iarg], "warn") != 0) { drtap = utils::numeric(FLERR, arg[iarg], false, lmp); iarg++; } } else if (strcmp(arg[iarg],"dsf") == 0) { vtype = 2; iarg++; - if (iarg < narg) { + // optional taper width; do not treat following keyword (e.g. "warn") as numeric + if (iarg < narg && strcmp(arg[iarg], "warn") != 0) { drtap = utils::numeric(FLERR, arg[iarg], false, lmp); iarg++; } From d100ea5c8d460396a1e649c359e19b3faa6ab454 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roman=20Gr=C3=B6ger?= <65243194+romangroger@users.noreply.github.com> Date: Tue, 27 Jan 2026 18:35:32 +0100 Subject: [PATCH 41/43] Update fix_qeq_slater.cpp --- src/QEQ/fix_qeq_slater.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/QEQ/fix_qeq_slater.cpp b/src/QEQ/fix_qeq_slater.cpp index c9dd0da19f7..35818ad05ee 100644 --- a/src/QEQ/fix_qeq_slater.cpp +++ b/src/QEQ/fix_qeq_slater.cpp @@ -58,7 +58,7 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg vtype = 1; iarg++; // optional taper width; do not treat following keyword (e.g. "warn") as numeric - if (iarg < narg && strcmp(arg[iarg], "warn") != 0) { + if (iarg < narg && utils::is_double(arg[iarg])) { drtap = utils::numeric(FLERR, arg[iarg], false, lmp); iarg++; } @@ -66,7 +66,7 @@ FixQEqSlater::FixQEqSlater(LAMMPS *lmp, int narg, char **arg) : FixQEq(lmp, narg vtype = 2; iarg++; // optional taper width; do not treat following keyword (e.g. "warn") as numeric - if (iarg < narg && strcmp(arg[iarg], "warn") != 0) { + if (iarg < narg && utils::is_double(arg[iarg])) { drtap = utils::numeric(FLERR, arg[iarg], false, lmp); iarg++; } From ae35daecbd90f670a76f30a57c4d6ea450321df7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 6 Feb 2026 14:49:37 -0500 Subject: [PATCH 42/43] small doc fixes and reformatting of doc sources --- doc/src/fix_qeq.rst | 22 +++++++++++++--- doc/src/pair_coul.rst | 28 +++++++++++++++++---- doc/utils/sphinx-config/false_positives.txt | 4 +-- 3 files changed, 43 insertions(+), 11 deletions(-) diff --git a/doc/src/fix_qeq.rst b/doc/src/fix_qeq.rst index ebc42d43885..464e159ea97 100644 --- a/doc/src/fix_qeq.rst +++ b/doc/src/fix_qeq.rst @@ -257,7 +257,21 @@ larger sizes, and *qeq/fire* is faster than *qeq/dynamic*\ . .. versionadded:: TBD -In previous versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. +In previous versions of LAMMPS, the real-space summations of Coulomb +interactions were done by replacing *1/r* using a damped potential +*erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of +decay. However, any finite value of *alpha* leads to a jump at the +cutoff, which interferes with equilibration if atoms move across the +cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` +(*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` +(*dsf*) solve this problem. An extension was implemented to specify the +width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the +Coulomb integrals at the cutoff. This is done by specifying the optional +arguments *wolf* and *dsf* with the value representing the width of +taper that smoothly terminates the Coulomb integrals. For example, if +the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are +smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For +backward compatibility, the default taper width is zero. Restart, fix_modify, output, run start/stop, minimize info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" @@ -325,14 +339,14 @@ Physical Chemistry, 105, 9396-9049 (2001) **(QEq/Fire)** T.-R. Shan, A. P. Thompson, S. J. Plimpton, in preparation -.. _Wolf: +.. _Wolf4: **(Wolf)** D. Wolf, P. Keblinski, S. R. Phillpot, J. Eggebrecht, J. Chem. Phys. 110, 8254 (1999). -.. _Fennell: +.. _Fennell3: **(Fennell)** J. Fennell, J. D. Gezelter, J. Chem. Phys. 124, 234104 (2006). -.. _Mei: +.. _Mei2: **(Mei)** J. Mei, J. W. Davenport, G. W. Fernando, Phys. Rev. B 43, 4653 (1991). diff --git a/doc/src/pair_coul.rst b/doc/src/pair_coul.rst index 205af881571..39cdaa4c740 100644 --- a/doc/src/pair_coul.rst +++ b/doc/src/pair_coul.rst @@ -308,9 +308,9 @@ command doc page. Alternatively *qfile* can be replaced by "coul/streitz", in which case the fix will extract QEq parameters from the coul/streitz pair style itself. -See the examples/streitz directory for an example input script that +See the ``examples/streitz`` directory for an example input script that uses the Streitz-Mintmire potential. The potentials directory has the -AlO.eam.alloy and AlO.streitz potential files used by the example. +``AlO.eam.alloy`` and ``AlO.streitz`` potential files used by the example. Note that the Streitz-Mintmire potential is generally used for oxides, but there is no conceptual problem with extending it to nitrides and @@ -321,9 +321,27 @@ Streitz-Mintmire parameterization for the material being modeled. .. versionchanged:: TBD -In previous versions of LAMMPS, the real-space summations of Coulomb interactions were done by replacing *1/r* using a damped potential *erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of decay. However, any finite value of *alpha* leads to a jump at the cutoff, which interferes with equilibration if atoms move across the cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` (*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` (*dsf*) solve this problem. An extension was implemented to specify the width of taper (see :ref:`(Mei et al.) `) to smoothly terminate the Coulomb integrals at the cutoff. This is done by specifying the optional arguments *wolf* and *dsf* with the value representing the width of taper that smoothly terminates the Coulomb integrals. For example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb integrals are smoothly rescaled from their actual value at r=6 A to zero at r=8 A. For backward compatibility, the default taper width is zero. - -An implementation of the Streitz-Mintmire potential for GaN due to :ref:`(Groger and Fikar) ` can be found in the examples/streitz directory. The electrostatic parameters of Ga and N are stored in file GaN.streitz and the short-range tersoff/mod potential in the file GaN.streitz+tersoff.mod. The total potential must be specified as: +In previous versions of LAMMPS, the real-space summations of Coulomb +interactions were done by replacing *1/r* using a damped potential +*erfc(alpha*r)/r* with the parameter *alpha* controlling the rate of +decay. However, any finite value of *alpha* leads to a jump at the +cutoff, which interferes with equilibration if atoms move across the +cutoff. The charge-neutralized potential of :ref:`(Wolf et al.) ` +(*wolf*) and its extension by :ref:`(Fennell and Gezelter) ` +(*dsf*) solve this problem. An extension was implemented to specify the +width of taper (see :ref:`(Mei et al.) `) to smoothly terminate +the Coulomb integrals at the cutoff. This is done by specifying the +optional arguments *wolf* and *dsf* with the value representing the +width of taper that smoothly terminates the Coulomb integrals. For +example, if the cutoff is 8 A and the taper width is 2 A, the Coulomb +integrals are smoothly rescaled from their actual value at r=6 A to zero +at r=8 A. For backward compatibility, the default taper width is zero. + +An implementation of the Streitz-Mintmire potential for GaN due to +:ref:`(Groger and Fikar) ` can be found in the examples/streitz +directory. The electrostatic parameters of Ga and N are stored in file +GaN.streitz and the short-range tersoff/mod potential in the file +GaN.streitz+tersoff.mod. The total potential must be specified as: .. code-block:: LAMMPS diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 8f2e1cf69d4..ffd35aec502 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1236,6 +1236,7 @@ Fickian fieldname figshare Fij +Fikar filelink filename Filename @@ -1464,6 +1465,7 @@ Grigera Grimme grmask Grmask +Groger gromacs Gromacs gromos @@ -3788,10 +3790,8 @@ Straub strcmp streitz Streitz -Streiz strerror strided -strietz stringstreams strmatch strncmp From 6032ea045e9187bc311915d4f36add140ce1a833 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 6 Feb 2026 14:51:53 -0500 Subject: [PATCH 43/43] replace potential files with symlinks --- examples/streitz/GaN.streitz | 4 +--- examples/streitz/GaN.streitz+tersoff.mod | 29 +----------------------- 2 files changed, 2 insertions(+), 31 deletions(-) mode change 100644 => 120000 examples/streitz/GaN.streitz mode change 100644 => 120000 examples/streitz/GaN.streitz+tersoff.mod diff --git a/examples/streitz/GaN.streitz b/examples/streitz/GaN.streitz deleted file mode 100644 index d5eb6cc929e..00000000000 --- a/examples/streitz/GaN.streitz +++ /dev/null @@ -1,3 +0,0 @@ -# chi (eV), J (eV), gamma (1/ang), zeta (1/ang), Z (e) -Ga -1.25 14.51 0 1.2333912356832242 0 -N 3.38 6.91 0 0.8509319105529678 0 diff --git a/examples/streitz/GaN.streitz b/examples/streitz/GaN.streitz new file mode 120000 index 00000000000..3a54d09561b --- /dev/null +++ b/examples/streitz/GaN.streitz @@ -0,0 +1 @@ +../../potentials/GaN.streitz \ No newline at end of file diff --git a/examples/streitz/GaN.streitz+tersoff.mod b/examples/streitz/GaN.streitz+tersoff.mod deleted file mode 100644 index 304ba73ceee..00000000000 --- a/examples/streitz/GaN.streitz+tersoff.mod +++ /dev/null @@ -1,28 +0,0 @@ -# DATE: 2026-01-27 CONTRIBUTOR: Roman Groger -# CITATION: R. Groger, J. Fikar, Acta Mater. XX, XXXX (2026) -# -# Modified Tersoff potential parameters for GaN describing short-range interactions in ionic-covalent -# model of GaN developed within the framework of the Streitz-Mintmire formalism. This file must be -# combined with the electrostatic energy to form a hybrid potential. The electrostatic parameters -# of Ga and N are stored in the file GaN.streitz. -# -# The parameters of homobonds are taken from Nord et al., J. Phys. Cond. Mat 15, 5649 (2003). -# Fitting of Ga-N parameters and the development of the ionic-covalent model for GaN is described in -# Groger & Fikar, Acta Mater. XX, XXXX (2026) -# -# Multiple entries can be added to this file, LAMMPS reads the ones it needs. -# These entries are in LAMMPS "metal" units: -# A,B = eV; lambda1,lambda2 = 1/Angstroms; R,D = Angstroms -# other quantities are unitless -# -# Format of a single entry (one or more lines): -# element 1, element 2, element 3, beta, alpha, h, eta, beta_ters, lambda2, B, R, D, lambda1, A, n, c1, c2, c3, c4, c5 - -Ga Ga Ga 1.0 1.846 -0.3013 1.0 1.0 1.4497 410.132 2.87 0.15 1.60916 535.199 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 -N N N 1.0 0.0 -0.045238 1.0 1.0 2.38426 423.769 2.2 0.2 3.55779 1044.77 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 -Ga Ga N 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -Ga N N 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -N Ga Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 2.7220181847973004 4626.922179188084 2.9 0.2 3.0086205179630774 7426.478709035391 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -N Ga N 1.0 0.0 -0.045238 1.0 1.0 0.0 0.0 2.2 0.2 0.0 0.0 1.0 0.76612 0.5998480604393894 0.040690958400000005 0.0 0.0 -N N Ga 1.0 0.0 -0.42139396832132586 1.0 1.0 0.0 0.0 2.9 0.2 0.0 0.0 1.0 0.0009358372029931801 0.07556968022478684 0.42407692241770484 0.0 0.0 -Ga N Ga 1.0 1.846 -0.3013 1.0 1.0 0.0 0.0 2.87 0.15 0.0 0.0 1.0 0.007874 0.05149559604622222 0.5625 0.0 0.0 diff --git a/examples/streitz/GaN.streitz+tersoff.mod b/examples/streitz/GaN.streitz+tersoff.mod new file mode 120000 index 00000000000..eb5a74aec23 --- /dev/null +++ b/examples/streitz/GaN.streitz+tersoff.mod @@ -0,0 +1 @@ +../../potentials/GaN.streitz+tersoff.mod \ No newline at end of file