diff --git a/src/GPU/pair_eam_gpu.cpp b/src/GPU/pair_eam_gpu.cpp index 1aac379c35..307d7795db 100644 --- a/src/GPU/pair_eam_gpu.cpp +++ b/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); } diff --git a/src/KOKKOS/pair_eam_kokkos.cpp b/src/KOKKOS/pair_eam_kokkos.cpp index 99a0756610..19f54796ba 100644 --- a/src/KOKKOS/pair_eam_kokkos.cpp +++ b/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::PairEAMKokkos(LAMMPS *lmp) : PairEAM(lmp) { respa_enable = 0; single_enable = 0; + unit_convert_flag = utils::NOCONVERT; atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 53f1b7a328..bd20d20ea7 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h index ed525a471c..6a8c9ed507 100644 --- a/src/MANYBODY/pair_eam.h +++ b/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 diff --git a/src/MANYBODY/pair_eam_alloy.cpp b/src/MANYBODY/pair_eam_alloy.cpp index 42587e9434..8f44492142 100644 --- a/src/MANYBODY/pair_eam_alloy.cpp +++ b/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(); diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp index 102e2a62f3..4eeeb173ec 100644 --- a/src/MANYBODY/pair_eam_cd.cpp +++ b/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); } } diff --git a/src/MANYBODY/pair_eam_fs.cpp b/src/MANYBODY/pair_eam_fs.cpp index 448a1d2701..4abf820861 100644 --- a/src/MANYBODY/pair_eam_fs.cpp +++ b/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(); diff --git a/src/USER-INTEL/pair_eam_alloy_intel.cpp b/src/USER-INTEL/pair_eam_alloy_intel.cpp index 16a55e4b2b..0f20b13440 100644 --- a/src/USER-INTEL/pair_eam_alloy_intel.cpp +++ b/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(); diff --git a/src/USER-INTEL/pair_eam_fs_intel.cpp b/src/USER-INTEL/pair_eam_fs_intel.cpp index e68d5687cf..5cc921e29e 100644 --- a/src/USER-INTEL/pair_eam_fs_intel.cpp +++ b/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(); diff --git a/src/USER-INTEL/pair_eam_intel.cpp b/src/USER-INTEL/pair_eam_intel.cpp index 994f0d3910..b1187c4433 100644 --- a/src/USER-INTEL/pair_eam_intel.cpp +++ b/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; diff --git a/src/USER-OMP/pair_eam_alloy_omp.cpp b/src/USER-OMP/pair_eam_alloy_omp.cpp index e8f56a9700..84f769a9ca 100644 --- a/src/USER-OMP/pair_eam_alloy_omp.cpp +++ b/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(); diff --git a/src/USER-OMP/pair_eam_fs_omp.cpp b/src/USER-OMP/pair_eam_fs_omp.cpp index c68d5f6019..96b4f77005 100644 --- a/src/USER-OMP/pair_eam_fs_omp.cpp +++ b/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(); diff --git a/src/USER-OMP/pair_eam_omp.cpp b/src/USER-OMP/pair_eam_omp.cpp index 60ae65def5..e5b0d0d7d2 100644 --- a/src/USER-OMP/pair_eam_omp.cpp +++ b/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); } diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp index abc02fedd4..f397846ffd 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/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