From e706cea897650b7f93d02b218be7bf50cd20311b Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 20 Sep 2012 16:03:32 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8807
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/USER-OMP/pair_eam_omp.cpp | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/USER-OMP/pair_eam_omp.cpp b/src/USER-OMP/pair_eam_omp.cpp
index 0d74a94ba0..59e58d1b73 100644
--- a/src/USER-OMP/pair_eam_omp.cpp
+++ b/src/USER-OMP/pair_eam_omp.cpp
@@ -184,6 +184,8 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
 
   // fp = derivative of embedding energy at each atom
   // phi = embedding energy at each atom
+  // if rho > rhomax (e.g. due to close approach of two atoms),
+  //   will exceed table, so add linear term to conserve energy
 
   for (ii = iifrom; ii < iito; ii++) {
     i = ilist[ii];
@@ -196,6 +198,7 @@ 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);
       e_tally_thr(this, i, i, nlocal, NEWTON_PAIR, phi, 0.0, thr);
     }
   }