From 3a94e4df2d2e981a1e4fd957947d16d8d5a85075 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 15 Apr 2024 02:50:16 -0400 Subject: [PATCH] print single warning when some rho[i] exceeded rhomax of the current EAM potential --- doc/src/pair_eam.rst | 14 ++++++++++++++ src/INTEL/pair_eam_intel.cpp | 18 +++++++++++++++--- src/MANYBODY/pair_eam.cpp | 18 +++++++++++++++++- src/MANYBODY/pair_eam.h | 3 ++- src/OPENMP/pair_eam_omp.cpp | 17 ++++++++++++++++- src/OPT/pair_eam_opt.cpp | 17 ++++++++++++++++- 6 files changed, 80 insertions(+), 7 deletions(-) diff --git a/doc/src/pair_eam.rst b/doc/src/pair_eam.rst index 654f49c166..bb6a01eabf 100644 --- a/doc/src/pair_eam.rst +++ b/doc/src/pair_eam.rst @@ -140,6 +140,20 @@ The OpenKIM Project at provides EAM potentials that can be used directly in LAMMPS with the :doc:`kim command ` interface. +.. warning:: + + The EAM potential files tabulate the embedding energy as a function + of the local electron density :math:`\rho`. When atoms get too + close, this electron density may exceed the range for which the + embedding energy was tabulated for. For simplicity and to avoid + errors during equilibration of randomized geometries, LAMMPS will + assume a linearly increasing embedding energy for electron densities + beyond the maximum tabulated value. This usually means that the EAM + model is not a good model for the kind of system under investigation. + LAMMPS will print a single warning when this happens. It may be + harmless at the beginning of an equilibration but would be a big + concern for accuracy if it happens during production runs. + ---------- For style *eam*, potential values are read from a file that is in the diff --git a/src/INTEL/pair_eam_intel.cpp b/src/INTEL/pair_eam_intel.cpp index bd78c3239d..4fd527858e 100644 --- a/src/INTEL/pair_eam_intel.cpp +++ b/src/INTEL/pair_eam_intel.cpp @@ -234,7 +234,7 @@ void PairEAMIntel::eval(const int offload, const int vflag, const int istride = fc.rhor_istride(); const int jstride = fc.rhor_jstride(); const int fstride = fc.frho_stride(); - + int beyond_rhomax = 0; { #if defined(__MIC__) && defined(_LMP_INTEL_OFFLOAD) *timer_compute = MIC_Wtime(); @@ -453,7 +453,10 @@ void PairEAMIntel::eval(const int offload, const int vflag, if (EFLAG) { flt_t phi = ((frho_spline_e[ioff].a*p + frho_spline_e[ioff].b)*p + frho_spline_e[ioff].c)*p + frho_spline_e[ioff].d; - if (rho[i] > frhomax) phi += fp_f[i] * (rho[i]-frhomax); + if (rho[i] > frhomax) { + phi += fp_f[i] * (rho[i]-frhomax); + beyond_rhomax = 1; + } if (!ONETYPE) { const int ptr_off=itype*ntypes + itype; oscale = scale_f[ptr_off]; @@ -653,6 +656,16 @@ void PairEAMIntel::eval(const int offload, const int vflag, else fix->stop_watch(TIME_HOST_PAIR); + if (EFLAG && (exceeded_rhomax >= 0)) { + MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world); + if (exceeded_rhomax > 0) { + if (comm->me == 0) + error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. " + "Computed embedding term is unreliable."); + exceeded_rhomax = -1; + } + } + if (EFLAG || vflag) fix->add_result_array(f_start, ev_global, offload, eatom, 0, vflag); else @@ -849,4 +862,3 @@ void PairEAMIntel::unpack_forward_comm(int n, int first, double *buf, last = first + n; for (i = first; i < last; i++) fp_f[i] = buf[m++]; } - diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 669a5cadbb..78690d38bf 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -43,6 +43,7 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp) unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nmax = 0; + exceeded_rhomax = 0; rho = nullptr; fp = nullptr; numforce = nullptr; @@ -145,6 +146,7 @@ void PairEAM::compute(int eflag, int vflag) double *coeff; int *ilist,*jlist,*numneigh,**firstneigh; + int beyond_rhomax = 0; evdwl = 0.0; ev_init(eflag,vflag); @@ -237,7 +239,9 @@ void PairEAM::compute(int eflag, int vflag) fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2]; if (eflag) { phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax); + if (rho[i] > rhomax) { + phi += fp[i] * (rho[i]-rhomax); + } phi *= scale[type[i]][type[i]]; if (eflag_global) eng_vdwl += phi; if (eflag_atom) eatom[i] += phi; @@ -322,6 +326,16 @@ void PairEAM::compute(int eflag, int vflag) } } + if (eflag && (exceeded_rhomax >= 0)) { + MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world); + if (exceeded_rhomax > 0) { + if (comm->me == 0) + error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. " + "Computed embedding term is unreliable."); + exceeded_rhomax = -1; + } + } + if (vflag_fdotr) virial_fdotr_compute(); } @@ -424,6 +438,8 @@ void PairEAM::init_style() neighbor->add_request(this); embedstep = -1; + + exceeded_rhomax = 0; } /* ---------------------------------------------------------------------- diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h index 67bbde570d..7ecc2ce0df 100644 --- a/src/MANYBODY/pair_eam.h +++ b/src/MANYBODY/pair_eam.h @@ -63,7 +63,8 @@ class PairEAM : public Pair { void swap_eam(double *, double **) override; protected: - int nmax; // allocated size of per-atom arrays + int nmax; // allocated size of per-atom arrays + int exceeded_rhomax; // 0 or 1 when rho[i] exceeded rhomax, -1 when not to check double cutforcesq; double **scale; bigint embedstep; // timestep, the embedding term was computed diff --git a/src/OPENMP/pair_eam_omp.cpp b/src/OPENMP/pair_eam_omp.cpp index e99fbedbb7..57aff9ea3a 100644 --- a/src/OPENMP/pair_eam_omp.cpp +++ b/src/OPENMP/pair_eam_omp.cpp @@ -17,6 +17,7 @@ #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" #include "neigh_list.h" @@ -102,6 +103,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr) double *coeff; int *ilist,*jlist,*numneigh,**firstneigh; + int beyond_rhomax = 0; evdwl = 0.0; const auto * _noalias const x = (dbl3_t *) atom->x[0]; @@ -203,7 +205,10 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr) fp[i] = (coeff[0]*p + coeff[1])*p + coeff[2]; if (EFLAG) { phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]; - if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax); + if (rho[i] > rhomax) { + phi += fp[i] * (rho[i]-rhomax); + beyond_rhomax = 1; + } e_tally_thr(this, i, i, nlocal, NEWTON_PAIR, scale[type[i]][type[i]]*phi, 0.0, thr); } } @@ -298,6 +303,16 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr) f[i].y += fytmp; f[i].z += fztmp; } + + if (EFLAG && (exceeded_rhomax >= 0)) { + MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world); + if (exceeded_rhomax > 0) { + if (comm->me == 0) + error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. " + "Computed embedding term is unreliable."); + exceeded_rhomax = -1; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/OPT/pair_eam_opt.cpp b/src/OPT/pair_eam_opt.cpp index 0560b0693a..6b4f1063e2 100644 --- a/src/OPT/pair_eam_opt.cpp +++ b/src/OPT/pair_eam_opt.cpp @@ -23,6 +23,7 @@ #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" #include "neigh_list.h" @@ -118,6 +119,7 @@ template void PairEAMOpt::eval() int ntypes = atom->ntypes; int ntypes2 = ntypes * ntypes; + int beyond_rhomax = 0; auto *_noalias fast_alpha = (fast_alpha_t *) malloc((size_t) ntypes2 * (nr + 1) * sizeof(fast_alpha_t)); @@ -251,7 +253,10 @@ template void PairEAMOpt::eval() fp[i] = (coeff[0] * p + coeff[1]) * p + coeff[2]; if (EFLAG) { double phi = ((coeff[3] * p + coeff[4]) * p + coeff[5]) * p + coeff[6]; - if (rho[i] > rhomax) phi += fp[i] * (rho[i] - rhomax); + if (rho[i] > rhomax) { + phi += fp[i] * (rho[i] - rhomax); + beyond_rhomax = 1; + } phi *= scale[type[i]][type[i]]; if (eflag_global) eng_vdwl += phi; if (eflag_atom) eatom[i] += phi; @@ -361,5 +366,15 @@ template void PairEAMOpt::eval() free(fast_gamma); fast_gamma = nullptr; + if (EFLAG && (exceeded_rhomax >= 0)) { + MPI_Allreduce(&beyond_rhomax, &exceeded_rhomax, 1, MPI_INT, MPI_MAX, world); + if (exceeded_rhomax > 0) { + if (comm->me == 0) + error->warning(FLERR, "Local rho[i] exceeded rhomax of EAM potential table. " + "Computed embedding term is unreliable."); + exceeded_rhomax = -1; + } + } + if (vflag_fdotr) virial_fdotr_compute(); }