lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

add unit conversion to pair styles eim, gw, gw/zbl, and nb3b/harmonic

This commit is contained in:
Axel Kohlmeyer 2020-06-28 16:27:32 -04:00
parent 9afbc71802
commit 435ac9078b
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
8 changed files with 241 additions and 33 deletions
Show Stats 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

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

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

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