add transparent unit conversion for most of the EAM family of potentials

src/GPU/pair_eam_gpu.cpp

@@ -31,6 +31,7 @@
 #include "neigh_request.h"
 #include "gpu_extra.h"
 #include "suffix.h"
+#include "utils.h"
 
 #define MAXLINE 1024
 
@@ -72,6 +73,7 @@ PairEAMGPU::PairEAMGPU(LAMMPS *lmp) : PairEAM(lmp), gpu_mode(GPU_FORCE)
   respa_enable = 0;
   reinitflag = 0;
   cpu_time = 0.0;
+  unit_convert_flag = utils::NOCONVERT;
   suffix_flag |= Suffix::GPU;
   GPU_EXTRA::gpu_ready(lmp->modify, lmp->error);
 }

src/KOKKOS/pair_eam_kokkos.cpp

@@ -28,6 +28,7 @@
 #include "memory_kokkos.h"
 #include "error.h"
 #include "atom_masks.h"
+#include "utils.h"
 
 using namespace LAMMPS_NS;
 
@@ -38,6 +39,7 @@ PairEAMKokkos<DeviceType>::PairEAMKokkos(LAMMPS *lmp) : PairEAM(lmp)
 {
   respa_enable = 0;
   single_enable = 0;
+  unit_convert_flag = utils::NOCONVERT;
 
   atomKK = (AtomKokkos *) atom;
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;

src/MANYBODY/pair_eam.cpp

@@ -43,6 +43,8 @@ PairEAM::PairEAM(LAMMPS *lmp) : Pair(lmp)
   restartinfo = 0;
   manybody_flag = 1;
   embedstep = -1;
+  unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
+  unit_convert_factor = 1.0;
 
   nmax = 0;
   rho = NULL;
@@ -241,7 +243,7 @@ void PairEAM::compute(int eflag, int vflag)
     if (eflag) {
       phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
       if (rho[i] > rhomax) phi += fp[i] * (rho[i]-rhomax);
-      phi *= scale[type[i]][type[i]];
+      phi *= scale[type[i]][type[i]] * unit_convert_factor;
       if (eflag_global) eng_vdwl += phi;
       if (eflag_atom) eatom[i] += phi;
     }
@@ -308,7 +310,7 @@ void PairEAM::compute(int eflag, int vflag)
         phi = z2*recip;
         phip = z2p*recip - phi*recip;
         psip = fp[i]*rhojp + fp[j]*rhoip + phip;
-        fpair = -scale[itype][jtype]*psip*recip;
+        fpair = -scale[itype][jtype]*psip*recip*unit_convert_factor;
 
         f[i][0] += delx*fpair;
         f[i][1] += dely*fpair;
@@ -319,7 +321,7 @@ void PairEAM::compute(int eflag, int vflag)
           f[j][2] -= delz*fpair;
         }
 
-        if (eflag) evdwl = scale[itype][jtype]*phi;
+        if (eflag) evdwl = scale[itype][jtype]*phi*unit_convert_factor;
         if (evflag) ev_tally(i,j,nlocal,newton_pair,
                              evdwl,0.0,fpair,delx,dely,delz);
       }
@@ -466,7 +468,12 @@ void PairEAM::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(lmp, filename, "EAM");
+    PotentialFileReader reader(lmp, filename, "EAM", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();
@@ -492,6 +499,7 @@ void PairEAM::read_file(char *filename)
       reader.next_dvector(&file->frho[1], file->nrho);
       reader.next_dvector(&file->zr[1], file->nr);
       reader.next_dvector(&file->rhor[1], file->nr);
+
     } catch (TokenizerException & e) {
       error->one(FLERR, e.what());
     }
@@ -834,9 +842,9 @@ double PairEAM::single(int i, int j, int itype, int jtype,
   phi += z2*recip;
   phip = z2p*recip - phi*recip;
   psip = fp[i]*rhojp + fp[j]*rhoip + phip;
-  fforce = -psip*recip;
+  fforce = -psip*recip*scale[i][j]*unit_convert_factor;
 
-  return phi;
+  return phi*scale[i][j]*unit_convert_factor;
 }
 
 /* ---------------------------------------------------------------------- */

src/MANYBODY/pair_eam.h

@@ -67,6 +67,7 @@ class PairEAM : public Pair {
   double cutforcesq;
   double **scale;
   bigint embedstep;           // timestep, the embedding term was computed
+  double unit_convert_factor; // multiplier for converting units transparently
 
   // per-atom arrays
 

src/MANYBODY/pair_eam_alloy.cpp

@@ -119,7 +119,12 @@ void PairEAMAlloy::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(lmp, filename, "EAM");
+    PotentialFileReader reader(lmp, filename, "EAMAlloy", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();

src/MANYBODY/pair_eam_cd.cpp

@@ -182,7 +182,7 @@ void PairEAMCD::compute(int eflag, int vflag)
     EAMTableIndex index = rhoToTableIndex(rho[i]);
     fp[i] = FPrimeOfRho(index, type[i]);
     if (eflag) {
-      phi = FofRho(index, type[i]);
+      phi = FofRho(index, type[i]) * unit_convert_factor;
       if (eflag_global) eng_vdwl += phi;
       if (eflag_atom) eatom[i] += phi;
     }
@@ -427,7 +427,7 @@ void PairEAMCD::compute(int eflag, int vflag)
 
         // Divide by r_ij and negate to get forces from gradient.
 
-        fpair /= -r;
+        fpair *= -unit_convert_factor/r;
 
         f[i][0] += delx*fpair;
         f[i][1] += dely*fpair;
@@ -438,7 +438,7 @@ void PairEAMCD::compute(int eflag, int vflag)
           f[j][2] -= delz*fpair;
         }
 
-        if (eflag) evdwl = phi;
+        if (eflag) evdwl = phi*unit_convert_factor;
         if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
       }
     }

src/MANYBODY/pair_eam_fs.cpp

@@ -119,7 +119,12 @@ void PairEAMFS::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(lmp, filename, "EAM");
+    PotentialFileReader reader(lmp, filename, "EAMFS", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();

src/USER-INTEL/pair_eam_alloy_intel.cpp

@@ -119,7 +119,12 @@ void PairEAMAlloyIntel::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(lmp, filename, "EAM");
+    PotentialFileReader reader(lmp, filename, "EAMAlloy", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();

src/USER-INTEL/pair_eam_fs_intel.cpp

@@ -119,7 +119,12 @@ void PairEAMFSIntel::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(lmp, filename, "EAM");
+    PotentialFileReader reader(lmp, filename, "EAMFS", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();

src/USER-INTEL/pair_eam_intel.cpp

@@ -435,7 +435,7 @@ void PairEAMIntel::eval(const int offload, const int vflag,
             const int ptr_off=itype*ntypes + itype;
             oscale = scale_f[ptr_off];
           }
-          phi *= oscale;
+          phi *= oscale * unit_convert_factor;
           tevdwl += phi;
           if (eatom) f[i].w += phi;
         }
@@ -549,7 +549,7 @@ void PairEAMIntel::eval(const int offload, const int vflag,
           const flt_t psip = fp_f[i]*rhojp + fp_f[j]*rhoip + phip;
           if (!ONETYPE)
             oscale = scale_fi[jtype];
-          const flt_t fpair = -oscale*psip*recip;
+          const flt_t fpair = -oscale*psip*recip*unit_convert_factor;
 
           const flt_t fpx = fpair * tdelx[jj];
           fxtmp += fpx;
@@ -562,7 +562,7 @@ void PairEAMIntel::eval(const int offload, const int vflag,
           if (NEWTON_PAIR) f[j].z -= fpz;
 
           if (EFLAG) {
-            const flt_t evdwl = oscale*phi;
+            const flt_t evdwl = oscale*phi*unit_convert_factor;
             sevdwl += evdwl;
             if (eatom) {
               fwtmp += (flt_t)0.5 * evdwl;

src/USER-OMP/pair_eam_alloy_omp.cpp

@@ -119,7 +119,13 @@ void PairEAMAlloyOMP::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(PairEAM::lmp, filename, "EAM");
+    PotentialFileReader reader(PairEAM::lmp, filename,
+                               "EAMAlloy", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();

src/USER-OMP/pair_eam_fs_omp.cpp

@@ -119,7 +119,13 @@ void PairEAMFSOMP::read_file(char *filename)
 
   // read potential file
   if(comm->me == 0) {
-    PotentialFileReader reader(PairEAM::lmp, filename, "EAM");
+    PotentialFileReader reader(PairEAM::lmp, filename,
+                               "EAMFS", unit_convert_flag);
+
+    // transparently convert units for supported conversions
+
+    unit_convert_factor
+      = utils::get_conversion_factor(utils::ENERGY, reader.get_unit_convert());
 
     try {
       reader.skip_line();

src/USER-OMP/pair_eam_omp.cpp

@@ -204,7 +204,8 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
     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, scale[type[i]][type[i]]*phi, 0.0, thr);
+      e_tally_thr(this, i, i, nlocal, NEWTON_PAIR,
+                  unit_convert_factor*scale[type[i]][type[i]]*phi, 0.0, thr);
     }
   }
 
@@ -278,7 +279,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
         phi = z2*recip;
         phip = z2p*recip - phi*recip;
         psip = fp[i]*rhojp + fp[j]*rhoip + phip;
-        fpair = -scale_i[jtype]*psip*recip;
+        fpair = -scale_i[jtype]*psip*recip*unit_convert_factor;
 
         fxtmp += delx*fpair;
         fytmp += dely*fpair;
@@ -289,7 +290,7 @@ void PairEAMOMP::eval(int iifrom, int iito, ThrData * const thr)
           f[j].z -= delz*fpair;
         }
 
-        if (EFLAG) evdwl = scale_i[jtype]*phi;
+        if (EFLAG) evdwl = scale_i[jtype]*phi*unit_convert_factor;
         if (EVFLAG) ev_tally_thr(this, i,j,nlocal,NEWTON_PAIR,
                                  evdwl,0.0,fpair,delx,dely,delz,thr);
       }

unittest/formats/test_pair_unit_convert.cpp

@@ -181,6 +181,49 @@ TEST_F(PairUnitConvertTest, lj_cut)
             EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
 }
 
+TEST_F(PairUnitConvertTest, eam)
+{
+    // check if the prerequisite pair style is available
+    if (!info->has_style("pair", "eam")) GTEST_SKIP();
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("units metal");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style eam");
+    lmp->input->one("pair_coeff * * Cu_u3.eam");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    // copy pressure, energy, and force from first step
+    double pold;
+    lmp->output->thermo->evaluate_keyword("press", &pold);
+    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    double **f  = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            fold[i][j] = f[i][j];
+
+    if (!verbose) ::testing::internal::CaptureStdout();
+    lmp->input->one("clear");
+    lmp->input->one("units real");
+    lmp->input->one("read_data test_pair_unit_convert.data");
+    lmp->input->one("pair_style eam");
+    lmp->input->one("pair_coeff * * Cu_u3.eam");
+    lmp->input->one("run 0 post no");
+    if (!verbose) ::testing::internal::GetCapturedStdout();
+
+    double pnew;
+    lmp->output->thermo->evaluate_keyword("press", &pnew);
+    EXPECT_NEAR(pold, p_convert * pnew, fabs(pnew * rel_error));
+    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
+    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));
+
+    f = lmp->atom->f;
+    for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 3; ++j)
+            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
+}
+
 TEST_F(PairUnitConvertTest, sw)
 {
     // check if the prerequisite pair style is available