diff --git a/doc/src/pair_eim.rst b/doc/src/pair_eim.rst index 1cea2f0993..107edda78d 100644 --- a/doc/src/pair_eim.rst +++ b/doc/src/pair_eim.rst @@ -48,7 +48,7 @@ and :math:`sigma_i` are calculated as \sigma_i = & \sum_{j=i_1}^{i_N} q_j \cdot \psi_{ij} \left(r_{ij}\right) \\ E_i\left(q_i,\sigma_i\right) = & \frac{1}{2} \cdot q_i \cdot \sigma_i -where :math:`\eta_{ji} is a pairwise function describing electron flow from atom +where :math:`\eta_{ji}` is a pairwise function describing electron flow from atom I to atom J, and :math:`\psi_{ij}` is another pairwise function. The multi-body nature of the EIM potential is a result of the embedding energy term. A complete list of all the pair functions used in EIM is summarized @@ -63,7 +63,7 @@ below \right.\\ \eta_{ji} = & A_{\eta,ij}\left(\chi_j-\chi_i\right)f_c\left(r,r_{s,\eta,ij},r_{c,\eta,ij}\right) \\ \psi_{ij}\left(r\right) = & A_{\psi,ij}\exp\left(-\zeta_{ij}r\right)f_c\left(r,r_{s,\psi,ij},r_{c,\psi,ij}\right) \\ - f_{c}\left(r,r_p,r_c\right) = & 0.510204 \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204 + f_{c}\left(r,r_p,r_c\right) = & 0.510204 \cdot \mathrm{erfc}\left[\frac{1.64498\left(2r-r_p-r_c\right)}{r_c-r_p}\right] - 0.010204 Here :math:`E_b, r_e, r_(c,\phi), \alpha, \beta, A_(\psi), \zeta, r_(s,\psi), r_(c,\psi), A_(\eta), r_(s,\eta), r_(c,\eta), \chi,` and pair function type diff --git a/examples/threebody/Si.sw b/examples/threebody/Si.sw deleted file mode 100644 index db4be100ef..0000000000 --- a/examples/threebody/Si.sw +++ /dev/null @@ -1,18 +0,0 @@ -# DATE: 2007-06-11 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Stillinger and Weber, Phys Rev B, 31, 5262, (1985) -# Stillinger-Weber parameters for various elements and mixtures -# multiple entries can be added to this file, LAMMPS reads the ones it needs -# these entries are in LAMMPS "metal" units: -# epsilon = eV; sigma = Angstroms -# other quantities are unitless - -# format of a single entry (one or more lines): -# element 1, element 2, element 3, -# epsilon, sigma, a, lambda, gamma, costheta0, A, B, p, q, tol - -# Here are the original parameters in metal units, for Silicon from: -# -# Stillinger and Weber, Phys. Rev. B, v. 31, p. 5262, (1985) -# - -Si Si Si 2.1683 2.0951 1.80 21.0 1.20 -0.333333333333 - 7.049556277 0.6022245584 4.0 0.0 0.0 diff --git a/examples/threebody/Si.sw b/examples/threebody/Si.sw new file mode 120000 index 0000000000..812023f52d --- /dev/null +++ b/examples/threebody/Si.sw @@ -0,0 +1 @@ +../../potentials/Si.sw \ No newline at end of file diff --git a/potentials/MOH.nb3b.harmonic b/potentials/MOH.nb3b.harmonic index 6bdc263b07..acf8f79c1c 100644 --- a/potentials/MOH.nb3b.harmonic +++ b/potentials/MOH.nb3b.harmonic @@ -1,4 +1,4 @@ -# DATE: 2016-07-29 UNITS: metal CONTRIBUTOR: Todd Zeitler, tzeitle@sandia.gov CITATION: none +# DATE: 2016-07-29 UNITS: real CONTRIBUTOR: Todd Zeitler, tzeitle@sandia.gov CITATION: none # nb3b/harmonic (nonbonded 3-body harmonic) parameters for various elements # # multiple entries can be added to this file, LAMMPS reads the ones it needs diff --git a/src/KOKKOS/sna_kokkos_impl.h b/src/KOKKOS/sna_kokkos_impl.h index 866b4c516f..d1de943fe3 100644 --- a/src/KOKKOS/sna_kokkos_impl.h +++ b/src/KOKKOS/sna_kokkos_impl.h @@ -552,10 +552,7 @@ void SNAKokkos::compute_bi(const int& iatom_mod, const int& jjb, con int idouble = 0; for (int elem1 = 0; elem1 < nelements; elem1++) { for (int elem2 = 0; elem2 < nelements; elem2++) { - const int jalloy = idouble; - for (int elem3 = 0; elem3 < nelements; elem3++) { - const int jjballoy = itriple; double sumzu = 0.0; double sumzu_temp = 0.0; @@ -566,7 +563,7 @@ void SNAKokkos::compute_bi(const int& iatom_mod, const int& jjb, con const int jjz_index = jjz+mb*(j+1)+ma; if (2*mb == j) return; // I think we can remove this? const auto utot = ulisttot_pack(iatom_mod, jju_index, elem3, iatom_div); - const auto zloc = zlist_pack(iatom_mod, jjz_index, jalloy, iatom_div); + const auto zloc = zlist_pack(iatom_mod, jjz_index, idouble, iatom_div); sumzu_temp += utot.re * zloc.re + utot.im * zloc.im; } } @@ -582,7 +579,7 @@ void SNAKokkos::compute_bi(const int& iatom_mod, const int& jjb, con const int jjz_index = jjz+(mb-1)*(j+1)+(j+1)+ma; const auto utot = ulisttot_pack(iatom_mod, jju_index, elem3, iatom_div); - const auto zloc = zlist_pack(iatom_mod, jjz_index, jalloy, iatom_div); + const auto zloc = zlist_pack(iatom_mod, jjz_index, idouble, iatom_div); sumzu_temp += utot.re * zloc.re + utot.im * zloc.im; } @@ -593,7 +590,7 @@ void SNAKokkos::compute_bi(const int& iatom_mod, const int& jjb, con const int jjz_index = jjz+(mb-1)*(j+1)+(j+1)+ma; const auto utot = ulisttot_pack(iatom_mod, jju_index, elem3, iatom_div); - const auto zloc = zlist_pack(iatom_mod, jjz_index, jalloy, iatom_div); + const auto zloc = zlist_pack(iatom_mod, jjz_index, idouble, iatom_div); sumzu += 0.5 * (utot.re * zloc.re + utot.im * zloc.im); } // end if jeven @@ -607,7 +604,7 @@ void SNAKokkos::compute_bi(const int& iatom_mod, const int& jjb, con sumzu -= bzero[j]; } } - blist_pack(iatom_mod, jjb, jjballoy, iatom_div) = sumzu; + blist_pack(iatom_mod, jjb, itriple, iatom_div) = sumzu; //} // end loop over j //} // end loop over j1, j2 itriple++; @@ -1078,7 +1075,7 @@ void SNAKokkos::compute_bi_cpu(const typename Kokkos::TeamPolicy::compute_bi_cpu(const typename Kokkos::TeamPolicyme == 0) { - EIMPotentialFileReader reader(lmp, filename); + EIMPotentialFileReader reader(lmp, filename, unit_convert_flag); reader.get_global(setfl); @@ -1050,14 +1051,18 @@ double PairEIM::memory_usage() return bytes; } -EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS * lmp, const std::string & filename) : +EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS *lmp, + const std::string &filename, + const int auto_convert) : Pointers(lmp), filename(filename) { if (comm->me != 0) { error->one(FLERR, "EIMPotentialFileReader should only be called by proc 0!"); } - FILE * fp = force->open_potential(filename.c_str()); + int unit_convert = auto_convert; + FILE *fp = force->open_potential(filename.c_str(), &unit_convert); + conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert); if (fp == NULL) { error->one(FLERR, fmt::format("cannot open EIM potential file {}", filename)); @@ -1186,7 +1191,7 @@ void EIMPotentialFileReader::parse(FILE * fp) PairData data; data.rcutphiA = values.next_double(); data.rcutphiR = values.next_double(); - data.Eb = values.next_double(); + data.Eb = values.next_double() * conversion_factor; data.r0 = values.next_double(); data.alpha = values.next_double(); data.beta = values.next_double(); @@ -1194,7 +1199,7 @@ void EIMPotentialFileReader::parse(FILE * fp) data.Asigma = values.next_double(); data.rq = values.next_double(); data.rcutsigma = values.next_double(); - data.Ac = values.next_double(); + data.Ac = values.next_double() * conversion_factor; data.zeta = values.next_double(); data.rs = values.next_double(); @@ -1217,20 +1222,18 @@ void EIMPotentialFileReader::parse(FILE * fp) } } -void EIMPotentialFileReader::get_global(PairEIM::Setfl * setfl) { +void EIMPotentialFileReader::get_global(PairEIM::Setfl *setfl) { setfl->division = division; setfl->rbig = rbig; setfl->rsmall = rsmall; } -void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const std::string & name) { - if (elements.find(name) == elements.end()) { - char str[128]; - snprintf(str, 128, "Element %s not defined in EIM potential file", name.c_str()); - error->one(FLERR, str); - } +void EIMPotentialFileReader::get_element(PairEIM::Setfl *setfl, int i, + const std::string &name) { + if (elements.find(name) == elements.end()) + error->one(FLERR,"Element " + name + " not defined in EIM potential file"); - ElementData & data = elements[name]; + ElementData &data = elements[name]; setfl->ielement[i] = data.ielement; setfl->mass[i] = data.mass; setfl->negativity[i] = data.negativity; @@ -1240,16 +1243,16 @@ void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const st setfl->q0[i] = data.q0; } -void EIMPotentialFileReader::get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB) { +void EIMPotentialFileReader::get_pair(PairEIM::Setfl *setfl, int ij, + const std::string &elemA, + const std::string &elemB) { auto p = get_pair(elemA, elemB); - if (pairs.find(p) == pairs.end()) { - char str[128]; - snprintf(str, 128, "Pair (%s, %s) not defined in EIM potential file", elemA.c_str(), elemB.c_str()); - error->one(FLERR, str); - } + if (pairs.find(p) == pairs.end()) + error->one(FLERR,"Element pair (" + elemA + ", " + elemB + + ") is not defined in EIM potential file"); - PairData & data = pairs[p]; + PairData &data = pairs[p]; setfl->rcutphiA[ij] = data.rcutphiA; setfl->rcutphiR[ij] = data.rcutphiR; setfl->Eb[ij] = data.Eb; diff --git a/src/MANYBODY/pair_eim.h b/src/MANYBODY/pair_eim.h index 986db52453..1266ffd242 100644 --- a/src/MANYBODY/pair_eim.h +++ b/src/MANYBODY/pair_eim.h @@ -98,17 +98,21 @@ class EIMPotentialFileReader : protected Pointers { std::string filename; static const int MAXLINE = 1024; char line[MAXLINE]; + double conversion_factor; - void parse(FILE * fp); - char * next_line(FILE * fp); - std::pair get_pair(const std::string & a, const std::string & b); + void parse(FILE *fp); + char *next_line(FILE *fp); + std::pair get_pair(const std::string &a, + const std::string &b); public: - EIMPotentialFileReader(class LAMMPS* lmp, const std::string & filename); + EIMPotentialFileReader(class LAMMPS* lmp, const std::string &filename, + const int auto_convert=0); - void get_global(PairEIM::Setfl * setfl); - void get_element(PairEIM::Setfl * setfl, int i, const std::string & name); - void get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB); + void get_global(PairEIM::Setfl *setfl); + void get_element(PairEIM::Setfl *setfl, int i, const std::string &name); + void get_pair(PairEIM::Setfl *setfl, int ij, + const std::string &elemA, const std::string &elemB); private: // potential parameters diff --git a/src/MANYBODY/pair_gw.cpp b/src/MANYBODY/pair_gw.cpp index e9604523b5..32b48f7eca 100644 --- a/src/MANYBODY/pair_gw.cpp +++ b/src/MANYBODY/pair_gw.cpp @@ -48,6 +48,7 @@ PairGW::PairGW(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; elements = NULL; @@ -375,9 +376,14 @@ void PairGW::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "GW"); + PotentialFileReader reader(lmp, file, "GW", unit_convert_flag); char * line; + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); while((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); @@ -427,6 +433,11 @@ void PairGW::read_file(char *file) params[nparams].lam1 = values.next_double(); params[nparams].biga = values.next_double(); params[nparams].powermint = int(params[nparams].powerm); + + if (unit_convert) { + params[nparams].biga *= conversion_factor; + params[nparams].bigb *= conversion_factor; + } } catch (TokenizerException & e) { error->one(FLERR, e.what()); } diff --git a/src/MANYBODY/pair_gw_zbl.cpp b/src/MANYBODY/pair_gw_zbl.cpp index 96fc742b42..ddc70174f6 100644 --- a/src/MANYBODY/pair_gw_zbl.cpp +++ b/src/MANYBODY/pair_gw_zbl.cpp @@ -70,9 +70,14 @@ void PairGWZBL::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "GW/ZBL"); + PotentialFileReader reader(lmp, file, "GW/ZBL", unit_convert_flag); char * line; + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); while((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); @@ -126,6 +131,11 @@ void PairGWZBL::read_file(char *file) params[nparams].ZBLcut = values.next_double(); params[nparams].ZBLexpscale = values.next_double(); params[nparams].powermint = int(params[nparams].powerm); + + if (unit_convert) { + params[nparams].biga *= conversion_factor; + params[nparams].bigb *= conversion_factor; + } } catch (TokenizerException & e) { error->one(FLERR, e.what()); } diff --git a/src/MANYBODY/pair_nb3b_harmonic.cpp b/src/MANYBODY/pair_nb3b_harmonic.cpp index 7491c375db..6f1bfb8905 100644 --- a/src/MANYBODY/pair_nb3b_harmonic.cpp +++ b/src/MANYBODY/pair_nb3b_harmonic.cpp @@ -47,6 +47,7 @@ PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; elements = NULL; @@ -291,9 +292,14 @@ void PairNb3bHarmonic::read_file(char *file) // open file on proc 0 if (comm->me == 0) { - PotentialFileReader reader(lmp, file, "nb3b/harmonic"); + PotentialFileReader reader(lmp, file, "nb3b/harmonic", unit_convert_flag); char * line; + // transparently convert units for supported conversions + + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY, + unit_convert); while((line = reader.next_line(NPARAMS_PER_LINE))) { try { ValueTokenizer values(line); @@ -331,6 +337,8 @@ void PairNb3bHarmonic::read_file(char *file) params[nparams].k_theta = values.next_double(); params[nparams].theta0 = values.next_double(); params[nparams].cutoff = values.next_double(); + + if (unit_convert) params[nparams].k_theta *= conversion_factor; } catch (TokenizerException & e) { error->one(FLERR, e.what()); } diff --git a/src/USER-REACTION/fix_bond_react.cpp b/src/USER-REACTION/fix_bond_react.cpp index 6ba5d1ce49..2c6c20c844 100644 --- a/src/USER-REACTION/fix_bond_react.cpp +++ b/src/USER-REACTION/fix_bond_react.cpp @@ -19,6 +19,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu) #include #include #include +#include #include "update.h" #include "modify.h" #include "respa.h" @@ -41,6 +42,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu) #include "error.h" #include "input.h" #include "variable.h" +#include "fmt/format.h" #include @@ -575,17 +577,11 @@ FixBondReact::~FixBondReact() delete [] set; if (group) { - char **newarg; - newarg = new char*[2]; - newarg[0] = master_group; - newarg[1] = (char *) "delete"; - group->assign(2,newarg); + group->assign(std::string(master_group) + " delete"); if (stabilization_flag == 1) { - newarg[0] = exclude_group; - group->assign(2,newarg); + group->assign(std::string(exclude_group) + " delete"); delete [] exclude_group; } - delete [] newarg; } } @@ -608,59 +604,38 @@ it will have the name 'i_limit_tags' and will be intitialized to 0 (not in group void FixBondReact::post_constructor() { + int len; // let's add the limit_tags per-atom property fix - int len = strlen("bond_react_props_internal") + 1; - id_fix2 = new char[len]; - strcpy(id_fix2,"bond_react_props_internal"); + std::string cmd = std::string("bond_react_props_internal"); + id_fix2 = new char[cmd.size()+1]; + strcpy(id_fix2,cmd.c_str()); int ifix = modify->find_fix(id_fix2); if (ifix == -1) { - char **newarg = new char*[7]; - newarg[0] = (char *) "bond_react_props_internal"; - newarg[1] = (char *) "all"; // group ID is ignored - newarg[2] = (char *) "property/atom"; - newarg[3] = (char *) "i_limit_tags"; - newarg[4] = (char *) "i_react_tags"; - newarg[5] = (char *) "ghost"; - newarg[6] = (char *) "yes"; - modify->add_fix(7,newarg); - delete [] newarg; + cmd += std::string(" all property/atom i_limit_tags i_react_tags ghost yes"); + modify->add_fix(cmd); } // create master_group if not already existing // NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart) group->find_or_create(master_group); - char **newarg; - newarg = new char*[5]; - newarg[0] = master_group; - newarg[1] = (char *) "dynamic"; - newarg[2] = (char *) "all"; - newarg[3] = (char *) "property"; - newarg[4] = (char *) "limit_tags"; - group->assign(5,newarg); - delete [] newarg; + cmd = fmt::format("{} dynamic all property limit_tags",master_group); + group->assign(cmd); if (stabilization_flag == 1) { int igroup = group->find(exclude_group); // create exclude_group if not already existing, or use as parent group if static if (igroup == -1 || group->dynamic[igroup] == 0) { // create stabilization per-atom property - len = strlen("bond_react_stabilization_internal") + 1; - id_fix3 = new char[len]; - strcpy(id_fix3,"bond_react_stabilization_internal"); + cmd = std::string("bond_react_stabilization_internal"); + id_fix3 = new char[cmd.size()+1]; + strcpy(id_fix3,cmd.c_str()); ifix = modify->find_fix(id_fix3); if (ifix == -1) { - char **newarg = new char*[6]; - newarg[0] = (char *) id_fix3; - newarg[1] = (char *) "all"; // group ID is ignored - newarg[2] = (char *) "property/atom"; - newarg[3] = (char *) "i_statted_tags"; - newarg[4] = (char *) "ghost"; - newarg[5] = (char *) "yes"; - modify->add_fix(6,newarg); + cmd += std::string(" all property/atom i_statted_tags ghost yes"); + modify->add_fix(cmd); fix3 = modify->fix[modify->nfix-1]; - delete [] newarg; } len = strlen("statted_tags") + 1; @@ -680,16 +655,11 @@ void FixBondReact::post_constructor() strcat(exclude_group,"_REACT"); group->find_or_create(exclude_group); - char **newarg; - newarg = new char*[5]; - newarg[0] = exclude_group; - newarg[1] = (char *) "dynamic"; - if (igroup == -1) newarg[2] = (char *) "all"; - else newarg[2] = (char *) exclude_PARENT_group; - newarg[3] = (char *) "property"; - newarg[4] = (char *) "statted_tags"; - group->assign(5,newarg); - delete [] newarg; + if (igroup == -1) + cmd = fmt::format("{} dynamic all property statted_tags",exclude_group); + else + cmd = fmt::format("{} dynamic {} property statted_tags",exclude_group,exclude_PARENT_group); + group->assign(cmd); delete [] exclude_PARENT_group; // on to statted_tags (system-wide thermostat) @@ -737,21 +707,16 @@ void FixBondReact::post_constructor() // let's create a new nve/limit fix to limit newly reacted atoms - len = strlen("bond_react_MASTER_nve_limit") + 1; - id_fix1 = new char[len]; - strcpy(id_fix1,"bond_react_MASTER_nve_limit"); + cmd = std::string("bond_react_MASTER_nve_limit"); + id_fix1 = new char[cmd.size()+1]; + strcpy(id_fix1,cmd.c_str()); ifix = modify->find_fix(id_fix1); if (ifix == -1) { - char **newarg = new char*[4]; - newarg[0] = id_fix1; - newarg[1] = master_group; - newarg[2] = (char *) "nve/limit"; - newarg[3] = nve_limit_xmax; - modify->add_fix(4,newarg); + cmd += fmt::format(" {} nve/limit {}",master_group,nve_limit_xmax); + modify->add_fix(cmd); fix1 = modify->fix[modify->nfix-1]; - delete [] newarg; } } } diff --git a/src/group.cpp b/src/group.cpp index be3ca5219a..61f28c2551 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -536,9 +536,9 @@ void Group::assign(int narg, char **arg) convenience function to allow assigning to groups from a single string ------------------------------------------------------------------------- */ -void Group::assign(const std::string &fixcmd) +void Group::assign(const std::string &groupcmd) { - std::vector args = utils::split_words(fixcmd); + std::vector args = utils::split_words(groupcmd); char **newarg = new char*[args.size()]; int i=0; for (const auto &arg : args) { diff --git a/src/pair_table.cpp b/src/pair_table.cpp index 73b916f966..d8f1d75958 100644 --- a/src/pair_table.cpp +++ b/src/pair_table.cpp @@ -1033,11 +1033,15 @@ void *PairTable::extract(const char *str, int &dim) if (strcmp(str,"cut_coul") != 0) return NULL; if (ntables == 0) error->all(FLERR,"All pair coeffs are not set"); - double cut_coul = tables[0].cut; - for (int m = 1; m < ntables; m++) - if (tables[m].cut != cut_coul) - error->all(FLERR, - "Pair table cutoffs must all be equal to use with KSpace"); - dim = 0; - return &tables[0].cut; + // only check for cutoff consistency if claiming to be KSpace compatible + + if (ewaldflag || pppmflag || msmflag || dispersionflag || tip4pflag) { + double cut_coul = tables[0].cut; + for (int m = 1; m < ntables; m++) + if (tables[m].cut != cut_coul) + error->all(FLERR, + "Pair table cutoffs must all be equal to use with KSpace"); + dim = 0; + return &tables[0].cut; + } else return NULL; } diff --git a/unittest/formats/test_pair_unit_convert.cpp b/unittest/formats/test_pair_unit_convert.cpp index b9e5431f6d..6c4b38397d 100644 --- a/unittest/formats/test_pair_unit_convert.cpp +++ b/unittest/formats/test_pair_unit_convert.cpp @@ -353,6 +353,178 @@ TEST_F(PairUnitConvertTest, eam_cd) EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error)); } +TEST_F(PairUnitConvertTest, eim) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "eim")) 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 eim"); + lmp->input->one("pair_coeff * * Na Cl ffield.eim Na Cl"); + 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 eim"); + lmp->input->one("pair_coeff * * Na Cl ffield.eim Na Cl"); + 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, gw) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "gw")) 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 gw"); + lmp->input->one("pair_coeff * * SiC.gw Si C"); + 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 gw"); + lmp->input->one("pair_coeff * * SiC.gw Si C"); + 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, gw_zbl) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "gw/zbl")) 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 gw/zbl"); + lmp->input->one("pair_coeff * * SiC.gw.zbl Si C"); + 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 gw/zbl"); + lmp->input->one("pair_coeff * * SiC.gw.zbl Si C"); + 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, nb3b_harmonic) +{ + // check if the prerequisite pair style is available + if (!info->has_style("pair", "nb3b/harmonic")) 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 nb3b/harmonic"); + lmp->input->one("pair_coeff * * MOH.nb3b.harmonic M O"); + 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 nb3b/harmonic"); + lmp->input->one("pair_coeff * * MOH.nb3b.harmonic M O"); + 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