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

add unit conversion to table pair style

Browse Source
This commit is contained in:
Axel Kohlmeyer
2020-06-25 04:31:34 -04:00
parent c6b16577f6
commit 2acf71c3e2
5 changed files with 140 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
src/pair_table.cpp
View File

@ -44,6 +44,7 @@ PairTable::PairTable(LAMMPS *lmp) : Pair(lmp)
{
ntables = 0;
tables = NULL;
unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);
}
/* ---------------------------------------------------------------------- */
@ -360,7 +361,13 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
{
TableFileReader reader(lmp, file, "pair");
char * line = reader.find_section_start(keyword);
// transparently convert units for supported conversions
int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,
unit_convert);
char *line = reader.find_section_start(keyword);
if (!line) {
error->one(FLERR,"Did not find keyword in table file");
@ -370,7 +377,7 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
// allocate table arrays for file values
line = reader.next_line();
param_extract(tb,line);
param_extract(tb, line);
memory->create(tb->rfile,tb->ninput,"pair:rfile");
memory->create(tb->efile,tb->ninput,"pair:efile");
memory->create(tb->ffile,tb->ninput,"pair:ffile");
@ -404,8 +411,8 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
ValueTokenizer values(line);
values.next_int();
rfile = values.next_double();
tb->efile[i] = values.next_double();
tb->ffile[i] = values.next_double();
tb->efile[i] = conversion_factor * values.next_double();
tb->ffile[i] = conversion_factor * values.next_double();
} catch (TokenizerException & e) {
++cerror;
}
@ -461,29 +468,26 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
}
}
if (ferror) {
std::string str = fmt::format("{} of {} force values in table are inconsistent with -dE/dr.\n"
" Should only be flagged at inflection points",ferror,tb->ninput);
error->warning(FLERR,str.c_str());
}
if (ferror)
error->warning(FLERR,fmt::format("{} of {} force values in table are "
"inconsistent with -dE/dr.\n Should "
"only be flagged at inflection points",
ferror,tb->ninput));
// warn if re-computed distance values differ from file values
if (rerror) {
char str[128];
sprintf(str,"%d of %d distance values in table with relative error\n"
" over %g to re-computed values",rerror,tb->ninput,EPSILONR);
error->warning(FLERR,str);
}
if (rerror)
error->warning(FLERR,fmt::format("{} of {} distance values in table with "
"relative error\n over {} to "
"re-computed values",
rerror,tb->ninput,EPSILONR));
// warn if data was read incompletely, e.g. columns were missing
if (cerror) {
char str[128];
sprintf(str,"%d of %d lines in table were incomplete\n"
" or could not be parsed completely",cerror,tb->ninput);
error->warning(FLERR,str);
}
if (cerror)
error->warning(FLERR,fmt::format("{} of {} lines in table were "
"incomplete\n or could not be parsed "
"completely",cerror,tb->ninput));
}
/* ----------------------------------------------------------------------

11
src/table_file_reader.cpp
View File

@ -27,17 +27,18 @@
using namespace LAMMPS_NS;
TableFileReader::TableFileReader(LAMMPS *lmp,
const std::string &filename,
const std::string &type) :
PotentialFileReader(lmp, filename, type + " table")
const std::string &filename,
const std::string &type,
const int auto_convert) :
PotentialFileReader(lmp, filename, type + " table", auto_convert)
{
}
TableFileReader::~TableFileReader() {
}
char * TableFileReader::find_section_start(const std::string & keyword) {
char * line = nullptr;
char *TableFileReader::find_section_start(const std::string & keyword) {
char *line = nullptr;
while ((line = reader->next_line())) {
ValueTokenizer values(line);
std::string word = values.next_string();

5
src/table_file_reader.h
View File

@ -24,10 +24,11 @@ namespace LAMMPS_NS
{
class TableFileReader : public PotentialFileReader {
public:
TableFileReader(class LAMMPS *lmp, const std::string &filename, const std::string & type);
TableFileReader(class LAMMPS *lmp, const std::string &filename,
const std::string &type, const int auto_convert = 0);
virtual ~TableFileReader();
char * find_section_start(const std::string & keyword);
char *find_section_start(const std::string &keyword);
};
} // namespace LAMMPS_NS

2
src/utils.cpp
View File

@ -687,7 +687,7 @@ double utils::get_conversion_factor(const int property, const int conversion)
} else if (conversion == METAL2REAL) {
return 23.060549;
} else if (conversion == REAL2METAL) {
return 0.04336410204284381954;
return 1.0/23.060549;
}
}
return 0.0;

105
unittest/formats/test_pair_unit_convert.cpp
View File

@ -140,6 +140,10 @@ TEST_F(PairUnitConvertTest, lj_cut)
lmp->input->one("read_data test_pair_unit_convert.data");
lmp->input->one("pair_style lj/cut 6.0");
lmp->input->one("pair_coeff * * 0.01014286346782117 2.0");
remove("test.table.metal");
lmp->input->one("pair_write 1 1 1000 r 0.1 6.0 test.table.metal lj_1_1");
lmp->input->one("pair_write 1 2 1000 r 0.1 6.0 test.table.metal lj_1_2");
lmp->input->one("pair_write 2 2 1000 r 0.1 6.0 test.table.metal lj_2_2");
lmp->input->one("run 0 post no");
if (!verbose) ::testing::internal::GetCapturedStdout();
@ -158,6 +162,10 @@ TEST_F(PairUnitConvertTest, lj_cut)
lmp->input->one("read_data test_pair_unit_convert.data");
lmp->input->one("pair_style lj/cut 6.0");
lmp->input->one("pair_coeff * * 0.2339 2.0");
remove("test.table.real");
lmp->input->one("pair_write 1 1 1000 r 0.1 6.0 test.table.real lj_1_1");
lmp->input->one("pair_write 1 2 1000 r 0.1 6.0 test.table.real lj_1_2");
lmp->input->one("pair_write 2 2 1000 r 0.1 6.0 test.table.real lj_2_2");
lmp->input->one("run 0 post no");
if (!verbose) ::testing::internal::GetCapturedStdout();
@ -216,6 +224,101 @@ TEST_F(PairUnitConvertTest, sw)
EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
}
TEST_F(PairUnitConvertTest, table_metal2real)
{
// check if the prerequisite pair style is available
if (!info->has_style("pair", "table")) 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 table linear 1000");
lmp->input->one("pair_coeff 1 1 test.table.metal lj_1_1");
lmp->input->one("pair_coeff 1 2 test.table.metal lj_1_2");
lmp->input->one("pair_coeff 2 2 test.table.metal lj_2_2");
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 table linear 1000");
lmp->input->one("pair_coeff 1 1 test.table.metal lj_1_1");
lmp->input->one("pair_coeff 1 2 test.table.metal lj_1_2");
lmp->input->one("pair_coeff 2 2 test.table.metal lj_2_2");
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, table_real2metal)
{
// check if the prerequisite pair style is available
if (!info->has_style("pair", "table")) GTEST_SKIP();
if (!verbose) ::testing::internal::CaptureStdout();
lmp->input->one("units real");
lmp->input->one("read_data test_pair_unit_convert.data");
lmp->input->one("pair_style table linear 1000");
lmp->input->one("pair_coeff 1 1 test.table.real lj_1_1");
lmp->input->one("pair_coeff 1 2 test.table.real lj_1_2");
lmp->input->one("pair_coeff 2 2 test.table.real lj_2_2");
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 metal");
lmp->input->one("read_data test_pair_unit_convert.data");
lmp->input->one("pair_style table linear 1000");
lmp->input->one("pair_coeff 1 1 test.table.real lj_1_1");
lmp->input->one("pair_coeff 1 2 test.table.real lj_1_2");
lmp->input->one("pair_coeff 2 2 test.table.real lj_2_2");
lmp->input->one("run 0 post no");
if (!verbose) ::testing::internal::GetCapturedStdout();
double pnew;
lmp->output->thermo->evaluate_keyword("press", &pnew);
EXPECT_NEAR(pold, 1.0/p_convert * pnew, fabs(pnew * rel_error));
double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
EXPECT_NEAR(1.0/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(1.0/ev_convert * fold[i][j], f[i][j],
fabs(f[i][j] * rel_error));
}
TEST_F(PairUnitConvertTest, tersoff)
{
// check if the prerequisite pair style is available
@ -539,5 +642,7 @@ int main(int argc, char **argv)
int rv = RUN_ALL_TESTS();
MPI_Finalize();
remove("test.table.metal");
remove("test.table.real");
return rv;
}