ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into more-internal-convenience

Browse Source 
# Conflicts:
#	src/group.cpp
This commit is contained in:
Axel Kohlmeyer
2020-06-29 11:50:43 -04:00
parent ab04416bfd d84a47707b
commit e6c4036a62
13 changed files with 288 additions and 131 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/pair_eim.rst
View File

@ -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

18
examples/threebody/Si.sw
View File

@ -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

1
examples/threebody/Si.sw Symbolic link
View File

@ -0,0 +1 @@
../../potentials/Si.sw

2
potentials/MOH.nb3b.harmonic
View File

@ -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

15
src/KOKKOS/sna_kokkos_impl.h
View File

@ -552,10 +552,7 @@ void SNAKokkos<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::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<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
int idouble = 0;
for (int elem1 = 0; elem1 < nelements; elem1++) {
for (int elem2 = 0; elem2 < nelements; elem2++) {
const auto jalloy = idouble;
auto jalloy = idouble; // must be non-const to work around gcc compiler bug
for (int elem3 = 0; elem3 < nelements; elem3++) {
Kokkos::parallel_for(Kokkos::TeamThreadRange(team,idxb_max),
[&] (const int& jjb) {
@ -1146,7 +1143,7 @@ void SNAKokkos<DeviceType>::compute_bi_cpu(const typename Kokkos::TeamPolicy<Dev
}
}
blist(jjb, jjballoy, iatom) = sumzu;
blist(jjb, itriple, iatom) = sumzu;
});
});
//} // end loop over j

43
src/MANYBODY/pair_eim.cpp
View File

@ -41,6 +41,7 @@ PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
setfl = NULL;
nmax = 0;
@ -477,7 +478,7 @@ void PairEIM::read_file(char *filename)
// read potential file
if( comm->me == 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;

18
src/MANYBODY/pair_eim.h
View File

@ -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<std::string, std::string> get_pair(const std::string & a, const std::string & b);
void parse(FILE *fp);
char *next_line(FILE *fp);
std::pair<std::string, std::string> 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

13
src/MANYBODY/pair_gw.cpp
View File

@ -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());
}

12
src/MANYBODY/pair_gw_zbl.cpp
View File

@ -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());
}

10
src/MANYBODY/pair_nb3b_harmonic.cpp
View File

@ -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());
}

89
src/USER-REACTION/fix_bond_react.cpp
View File

@ -19,6 +19,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
#include <mpi.h>
#include <cmath>
#include <cstring>
#include <string>
#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 <algorithm>
@ -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;
}
}
}

4
src/group.cpp
View File

@ -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<std::string> args = utils::split_words(fixcmd);
std::vector<std::string> args = utils::split_words(groupcmd);
char **newarg = new char*[args.size()];
int i=0;
for (const auto &arg : args) {

18
src/pair_table.cpp
View File

@ -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;
}

172
unittest/formats/test_pair_unit_convert.cpp
View File

@ -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