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

Update pair gw zbl

Browse Source
This commit is contained in:
Richard Berger
2020-05-28 15:56:36 -04:00
parent 607581e7f3
commit e206647717
3 changed files with 92 additions and 122 deletions
Download Patch File Download Diff File Expand all files Collapse all files

201
src/MANYBODY/pair_gw_zbl.cpp
View File

@ -28,6 +28,8 @@
#include "memory.h" #include "memory.h"
#include "error.h" #include "error.h"
#include "utils.h" #include "utils.h"
#include "tokenizer.h"
#include "potential_file_reader.h"
#include "math_const.h" #include "math_const.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -61,146 +63,101 @@ PairGWZBL::PairGWZBL(LAMMPS *lmp) : PairGW(lmp)
void PairGWZBL::read_file(char *file) void PairGWZBL::read_file(char *file)
{ {
int params_per_line = 21;
char **words = new char*[params_per_line+1];
memory->sfree(params); memory->sfree(params);
params = NULL; params = NULL;
nparams = maxparam = 0; nparams = maxparam = 0;
// open file on proc 0 // open file on proc 0
FILE *fp;
if (comm->me == 0) { if (comm->me == 0) {
fp = force->open_potential(file); PotentialFileReader reader(lmp, file, "GW/ZBL");
if (fp == NULL) { char * line;
char str[128];
snprintf(str,128,"Cannot open GW potential file %s",file);
error->one(FLERR,str);
}
}
// read each line out of file, skipping blank lines or leading '#' while(line = reader.next_line(NPARAMS_PER_LINE)) {
// store line of params if all 3 element tags are in element list try {
ValueTokenizer values(line, " \t\n\r\f");
int n,nwords,ielement,jelement,kelement; std::string iname = values.next_string();
char line[MAXLINE],*ptr; std::string jname = values.next_string();
int eof = 0; std::string kname = values.next_string();
while (1) { // ielement,jelement,kelement = 1st args
if (comm->me == 0) { // if all 3 args are in element list, then parse this line
ptr = fgets(line,MAXLINE,fp); // else skip to next entry in file
if (ptr == NULL) { int ielement, jelement, kelement;
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank for (ielement = 0; ielement < nelements; ielement++)
if (iname == elements[ielement]) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (jname == elements[jelement]) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (kname == elements[kelement]) break;
if (kelement == nelements) continue;
if ((ptr = strchr(line,'#'))) *ptr = '\0'; // load up parameter settings and error check their values
nwords = utils::count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
while (nwords < params_per_line) { params[nparams].ielement = ielement;
n = strlen(line); params[nparams].jelement = jelement;
if (comm->me == 0) { params[nparams].kelement = kelement;
ptr = fgets(&line[n],MAXLINE-n,fp); params[nparams].powerm = values.next_double();
if (ptr == NULL) { params[nparams].gamma = values.next_double();
eof = 1; params[nparams].lam3 = values.next_double();
fclose(fp); params[nparams].c = values.next_double();
} else n = strlen(line) + 1; params[nparams].d = values.next_double();
params[nparams].h = values.next_double();
params[nparams].powern = values.next_double();
params[nparams].beta = values.next_double();
params[nparams].lam2 = values.next_double();
params[nparams].bigb = values.next_double();
params[nparams].bigr = values.next_double();
params[nparams].bigd = values.next_double();
params[nparams].lam1 = values.next_double();
params[nparams].biga = values.next_double();
params[nparams].Z_i = values.next_double();
params[nparams].Z_j = values.next_double();
params[nparams].ZBLcut = values.next_double();
params[nparams].ZBLexpscale = values.next_double();
params[nparams].powermint = int(params[nparams].powerm);
} catch (TokenizerException & e) {
error->one(FLERR, e.what());
} }
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break; // currently only allow m exponent of 1 or 3
MPI_Bcast(&n,1,MPI_INT,0,world); if (
MPI_Bcast(line,n,MPI_CHAR,0,world); params[nparams].lam3 < 0.0 || params[nparams].c < 0.0 ||
if ((ptr = strchr(line,'#'))) *ptr = '\0'; params[nparams].d < 0.0 || params[nparams].powern < 0.0 ||
nwords = utils::count_words(line); params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||
params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
params[nparams].bigd < 0.0 ||
params[nparams].bigd > params[nparams].bigr ||
params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].powerm - params[nparams].powermint != 0.0 ||
(params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
params[nparams].gamma < 0.0 ||
params[nparams].Z_i < 1.0 || params[nparams].Z_j < 1.0 ||
params[nparams].ZBLcut < 0.0 || params[nparams].ZBLexpscale < 0.0)
error->one(FLERR,"Illegal GW parameter");
nparams++;
} }
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in GW potential file");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next line
for (ielement = 0; ielement < nelements; ielement++)
if (strcmp(words[0],elements[ielement]) == 0) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (strcmp(words[1],elements[jelement]) == 0) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (strcmp(words[2],elements[kelement]) == 0) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].powerm = atof(words[3]);
params[nparams].gamma = atof(words[4]);
params[nparams].lam3 = atof(words[5]);
params[nparams].c = atof(words[6]);
params[nparams].d = atof(words[7]);
params[nparams].h = atof(words[8]);
params[nparams].powern = atof(words[9]);
params[nparams].beta = atof(words[10]);
params[nparams].lam2 = atof(words[11]);
params[nparams].bigb = atof(words[12]);
params[nparams].bigr = atof(words[13]);
params[nparams].bigd = atof(words[14]);
params[nparams].lam1 = atof(words[15]);
params[nparams].biga = atof(words[16]);
params[nparams].Z_i = atof(words[17]);
params[nparams].Z_j = atof(words[18]);
params[nparams].ZBLcut = atof(words[19]);
params[nparams].ZBLexpscale = atof(words[20]);
// currently only allow m exponent of 1 or 3
params[nparams].powermint = int(params[nparams].powerm);
if (
params[nparams].lam3 < 0.0 || params[nparams].c < 0.0 ||
params[nparams].d < 0.0 || params[nparams].powern < 0.0 ||
params[nparams].beta < 0.0 || params[nparams].lam2 < 0.0 ||
params[nparams].bigb < 0.0 || params[nparams].bigr < 0.0 ||
params[nparams].bigd < 0.0 ||
params[nparams].bigd > params[nparams].bigr ||
params[nparams].lam3 < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].powerm - params[nparams].powermint != 0.0 ||
(params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
params[nparams].gamma < 0.0 ||
params[nparams].Z_i < 1.0 || params[nparams].Z_j < 1.0 ||
params[nparams].ZBLcut < 0.0 || params[nparams].ZBLexpscale < 0.0)
error->all(FLERR,"Illegal GW parameter");
nparams++;
} }
delete [] words; MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if(comm->me != 0) {
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

2
src/MANYBODY/pair_gw_zbl.h
View File

@ -29,6 +29,8 @@ class PairGWZBL : public PairGW {
PairGWZBL(class LAMMPS *); PairGWZBL(class LAMMPS *);
~PairGWZBL() {} ~PairGWZBL() {}
static const int NPARAMS_PER_LINE = 21;
private: private:
double global_a_0; // Bohr radius for Coulomb repulsion double global_a_0; // Bohr radius for Coulomb repulsion
double global_epsilon_0; // permittivity of vacuum for Coulomb repulsion double global_epsilon_0; // permittivity of vacuum for Coulomb repulsion

11
unittest/formats/test_potential_file_reader.cpp
View File

@ -11,6 +11,7 @@
#include "MANYBODY/pair_tersoff_mod_c.h" #include "MANYBODY/pair_tersoff_mod_c.h"
#include "MANYBODY/pair_tersoff_zbl.h" #include "MANYBODY/pair_tersoff_zbl.h"
#include "MANYBODY/pair_gw.h" #include "MANYBODY/pair_gw.h"
#include "MANYBODY/pair_gw_zbl.h"
#include <mpi.h> #include <mpi.h>
@ -24,6 +25,7 @@ const int LAMMPS_NS::PairTersoffMOD::NPARAMS_PER_LINE;
const int LAMMPS_NS::PairTersoffMODC::NPARAMS_PER_LINE; const int LAMMPS_NS::PairTersoffMODC::NPARAMS_PER_LINE;
const int LAMMPS_NS::PairTersoffZBL::NPARAMS_PER_LINE; const int LAMMPS_NS::PairTersoffZBL::NPARAMS_PER_LINE;
const int LAMMPS_NS::PairGW::NPARAMS_PER_LINE; const int LAMMPS_NS::PairGW::NPARAMS_PER_LINE;
const int LAMMPS_NS::PairGWZBL::NPARAMS_PER_LINE;
class PotenialFileReaderTest : public ::testing::Test { class PotenialFileReaderTest : public ::testing::Test {
protected: protected:
@ -117,6 +119,15 @@ TEST_F(PotenialFileReaderTest, GW) {
ASSERT_EQ(utils::count_words(line), PairGW::NPARAMS_PER_LINE); ASSERT_EQ(utils::count_words(line), PairGW::NPARAMS_PER_LINE);
} }
TEST_F(PotenialFileReaderTest, GWZBL) {
::testing::internal::CaptureStdout();
PotentialFileReader reader(lmp, "SiC.gw.zbl", "GWZBL");
::testing::internal::GetCapturedStdout();
auto line = reader.next_line(PairGWZBL::NPARAMS_PER_LINE);
ASSERT_EQ(utils::count_words(line), PairGWZBL::NPARAMS_PER_LINE);
}
int main(int argc, char **argv) int main(int argc, char **argv)
{ {
MPI_Init(&argc, &argv); MPI_Init(&argc, &argv);