Update pair eim
This commit is contained in:
@ -27,11 +27,11 @@
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "utils.h"
|
||||
#include "tokenizer.h"
|
||||
#include "potential_file_reader.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MAXLINE 1024
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
|
||||
@ -451,20 +451,6 @@ double PairEIM::init_one(int i, int j)
|
||||
|
||||
void PairEIM::read_file(char *filename)
|
||||
{
|
||||
// open potential file
|
||||
|
||||
int me = comm->me;
|
||||
FILE *fptr;
|
||||
|
||||
if (me == 0) {
|
||||
fptr = force->open_potential(filename);
|
||||
if (fptr == NULL) {
|
||||
char str[128];
|
||||
snprintf(str,128,"Cannot open EIM potential file %s",filename);
|
||||
error->one(FLERR,str);
|
||||
}
|
||||
}
|
||||
|
||||
int npair = nelements*(nelements+1)/2;
|
||||
setfl->ielement = new int[nelements];
|
||||
setfl->mass = new double[nelements];
|
||||
@ -488,52 +474,55 @@ void PairEIM::read_file(char *filename)
|
||||
setfl->rs = new double[npair];
|
||||
setfl->tp = new int[npair];
|
||||
|
||||
if (me == 0)
|
||||
if (!grabglobal(fptr))
|
||||
error->one(FLERR,"Could not grab global entry from EIM potential file");
|
||||
MPI_Bcast(&setfl->division,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rbig,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rsmall,1,MPI_DOUBLE,0,world);
|
||||
// read potential file
|
||||
if( comm->me == 0) {
|
||||
EIMPotentialFileReader reader(lmp, filename);
|
||||
|
||||
for (int i = 0; i < nelements; i++) {
|
||||
if (me == 0)
|
||||
if (!grabsingle(fptr,i))
|
||||
error->one(FLERR,"Could not grab element entry from EIM potential file");
|
||||
MPI_Bcast(&setfl->ielement[i],1,MPI_INT,0,world);
|
||||
MPI_Bcast(&setfl->mass[i],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->negativity[i],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->ra[i],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->ri[i],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->Ec[i],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->q0[i],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
reader.get_global(setfl);
|
||||
|
||||
for (int i = 0; i < nelements; i++) {
|
||||
for (int j = i; j < nelements; j++) {
|
||||
int ij;
|
||||
if (i == j) ij = i;
|
||||
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
|
||||
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
|
||||
if (me == 0)
|
||||
if (grabpair(fptr,i,j) == 0)
|
||||
error->one(FLERR,"Could not grab pair entry from EIM potential file");
|
||||
MPI_Bcast(&setfl->rcutphiA[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rcutphiR[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->Eb[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->r0[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->alpha[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->beta[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rcutq[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->Asigma[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rq[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rcutsigma[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->Ac[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->zeta[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->rs[ij],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&setfl->tp[ij],1,MPI_INT,0,world);
|
||||
for (int i = 0; i < nelements; i++) {
|
||||
reader.get_element(setfl, i, elements[i]);
|
||||
}
|
||||
|
||||
for (int i = 0; i < nelements; i++) {
|
||||
for (int j = i; j < nelements; j++) {
|
||||
int ij;
|
||||
if (i == j) ij = i;
|
||||
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
|
||||
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
|
||||
reader.get_pair(setfl, ij, elements[i], elements[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// broadcast potential information to other procs
|
||||
MPI_Bcast(&setfl->division, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&setfl->rbig, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&setfl->rsmall, 1, MPI_DOUBLE, 0, world);
|
||||
|
||||
MPI_Bcast(setfl->ielement, nelements, MPI_INT, 0, world);
|
||||
MPI_Bcast(setfl->mass, nelements, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->negativity, nelements, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->ra, nelements, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->ri, nelements, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->Ec, nelements, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->q0, nelements, MPI_DOUBLE, 0, world);
|
||||
|
||||
MPI_Bcast(setfl->rcutphiA, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->rcutphiR, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->Eb, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->r0, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->alpha, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->beta, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->rcutq, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->Asigma, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->rq, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->rcutsigma, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->Ac, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->zeta, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->rs, npair, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(setfl->tp, npair, MPI_INT, 0, world);
|
||||
|
||||
setfl->nr = 5000;
|
||||
setfl->cut = 0.0;
|
||||
for (int i = 0; i < npair; i++) {
|
||||
@ -602,10 +591,6 @@ void PairEIM::read_file(char *filename)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// close the potential file
|
||||
|
||||
if (me == 0) fclose(fptr);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -879,114 +864,6 @@ void PairEIM::interpolate(int n, double delta, double *f,
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
grab global line from file and store info in setfl
|
||||
return 0 if error
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int PairEIM::grabglobal(FILE *fptr)
|
||||
{
|
||||
char line[MAXLINE];
|
||||
|
||||
char *pch = NULL, *data = NULL;
|
||||
while (pch == NULL) {
|
||||
if (fgets(line,MAXLINE,fptr) == NULL) break;
|
||||
pch = strstr(line,"global");
|
||||
if (pch != NULL) {
|
||||
data = strtok (line," \t\n\r\f");
|
||||
data = strtok (NULL,"?");
|
||||
sscanf(data,"%lg %lg %lg",&setfl->division,&setfl->rbig,&setfl->rsmall);
|
||||
}
|
||||
}
|
||||
if (pch == NULL) return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
grab elemental line from file and store info in setfl
|
||||
return 0 if error
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int PairEIM::grabsingle(FILE *fptr, int i)
|
||||
{
|
||||
char line[MAXLINE];
|
||||
|
||||
rewind(fptr);
|
||||
|
||||
char *pch1 = NULL, *pch2 = NULL, *data = NULL;
|
||||
while (pch1 == NULL || pch2 == NULL) {
|
||||
if (fgets(line,MAXLINE,fptr) == NULL) break;
|
||||
pch1 = strtok (line," \t\n\r\f");
|
||||
pch1 = strstr(pch1,"element:");
|
||||
if (pch1 != NULL) {
|
||||
pch2 = strtok(NULL, " \t\n\r\f");
|
||||
if (pch2 != NULL) {
|
||||
data = strtok (NULL, "?");
|
||||
if (strcmp(pch2,elements[i]) == 0) {
|
||||
sscanf(data,"%d %lg %lg %lg %lg %lg %lg",&setfl->ielement[i],
|
||||
&setfl->mass[i],&setfl->negativity[i],&setfl->ra[i],
|
||||
&setfl->ri[i],&setfl->Ec[i],&setfl->q0[i]);
|
||||
} else pch2 = NULL;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (pch1 == NULL || pch2 == NULL) return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
grab pair line from file and store info in setfl
|
||||
return 0 if error
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
int PairEIM::grabpair(FILE *fptr, int i, int j)
|
||||
{
|
||||
char line[MAXLINE];
|
||||
|
||||
rewind(fptr);
|
||||
|
||||
int ij;
|
||||
if (i == j) ij = i;
|
||||
else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
|
||||
else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
|
||||
|
||||
char *pch1 = NULL, *pch2 = NULL, *pch3 = NULL, *data = NULL;
|
||||
while (pch1 == NULL || pch2 == NULL || pch3 == NULL) {
|
||||
if (fgets(line,MAXLINE,fptr) == NULL) break;
|
||||
pch1 = strtok (line," \t\n\r\f");
|
||||
pch1 = strstr(pch1,"pair:");
|
||||
if (pch1 != NULL) {
|
||||
pch2 = strtok (NULL, " \t\n\r\f");
|
||||
if (pch2 != NULL) pch3 = strtok (NULL, " \t\n\r\f");
|
||||
if (pch3 != NULL) data = strtok (NULL, "?");
|
||||
if ((pch2 != NULL) && (pch3 != NULL)) {
|
||||
if ((strcmp(pch2,elements[i]) == 0 &&
|
||||
strcmp(pch3,elements[j]) == 0) ||
|
||||
(strcmp(pch2,elements[j]) == 0 &&
|
||||
strcmp(pch3,elements[i]) == 0)) {
|
||||
sscanf(data,"%lg %lg %lg %lg %lg",
|
||||
&setfl->rcutphiA[ij],&setfl->rcutphiR[ij],
|
||||
&setfl->Eb[ij],&setfl->r0[ij],&setfl->alpha[ij]);
|
||||
utils::sfgets(FLERR,line,MAXLINE,fptr,NULL,error);
|
||||
sscanf(line,"%lg %lg %lg %lg %lg",
|
||||
&setfl->beta[ij],&setfl->rcutq[ij],&setfl->Asigma[ij],
|
||||
&setfl->rq[ij],&setfl->rcutsigma[ij]);
|
||||
utils::sfgets(FLERR,line,MAXLINE,fptr,NULL,error);
|
||||
sscanf(line,"%lg %lg %lg %d",
|
||||
&setfl->Ac[ij],&setfl->zeta[ij],&setfl->rs[ij],
|
||||
&setfl->tp[ij]);
|
||||
} else {
|
||||
pch1 = NULL;
|
||||
pch2 = NULL;
|
||||
pch3 = NULL;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (pch1 == NULL || pch2 == NULL || pch3 == NULL) return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
cutoff function
|
||||
------------------------------------------------------------------------- */
|
||||
@ -1171,3 +1048,221 @@ double PairEIM::memory_usage()
|
||||
bytes += 2 * nmax * sizeof(double);
|
||||
return bytes;
|
||||
}
|
||||
|
||||
EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS * lmp, const std::string & filename) :
|
||||
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());
|
||||
|
||||
if (fp == NULL) {
|
||||
char str[128];
|
||||
snprintf(str, 128, "cannot open EIM potential file %s", filename.c_str());
|
||||
error->one(FLERR, str);
|
||||
}
|
||||
|
||||
parse(fp);
|
||||
|
||||
fclose(fp);
|
||||
}
|
||||
|
||||
std::pair<std::string, std::string> EIMPotentialFileReader::get_pair(const std::string & a, const std::string & b) {
|
||||
if (a < b) {
|
||||
return std::make_pair(a, b);
|
||||
}
|
||||
return std::make_pair(b, a);
|
||||
}
|
||||
|
||||
char * EIMPotentialFileReader::next_line(FILE * fp) {
|
||||
// concatenate lines if they end with '&'
|
||||
// strip comments after '#'
|
||||
int n = 0;
|
||||
int nwords = 0;
|
||||
bool concat = false;
|
||||
|
||||
char *ptr = fgets(line, MAXLINE, fp);
|
||||
|
||||
if (ptr == nullptr) {
|
||||
// EOF
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
// strip comment
|
||||
if ((ptr = strchr(line, '#'))) *ptr = '\0';
|
||||
|
||||
// strip ampersand
|
||||
if ((ptr = strrchr(line, '&'))) {
|
||||
concat = true;
|
||||
*ptr = '\0';
|
||||
}
|
||||
|
||||
nwords = utils::count_words(line);
|
||||
|
||||
if (nwords > 0) {
|
||||
n = strlen(line);
|
||||
}
|
||||
|
||||
while(n == 0 || concat) {
|
||||
char *ptr = fgets(&line[n], MAXLINE - n, fp);
|
||||
|
||||
if (ptr == nullptr) {
|
||||
// EOF
|
||||
return line;
|
||||
}
|
||||
|
||||
// strip comment
|
||||
if ((ptr = strchr(line, '#'))) *ptr = '\0';
|
||||
|
||||
// strip ampersand
|
||||
if ((ptr = strrchr(line, '&'))) {
|
||||
concat = true;
|
||||
*ptr = '\0';
|
||||
} else {
|
||||
concat = false;
|
||||
}
|
||||
|
||||
nwords = utils::count_words(line);
|
||||
|
||||
// skip line if blank
|
||||
if (nwords > 0) {
|
||||
n = strlen(line);
|
||||
}
|
||||
}
|
||||
|
||||
return line;
|
||||
}
|
||||
|
||||
void EIMPotentialFileReader::parse(FILE * fp)
|
||||
{
|
||||
char * line = nullptr;
|
||||
bool found_global = false;
|
||||
|
||||
while(line = next_line(fp)) {
|
||||
ValueTokenizer values(line, " \t\r\n\f");
|
||||
std::string type = values.next_string();
|
||||
|
||||
if (type == "global:") {
|
||||
if (values.count() != 4) {
|
||||
error->one(FLERR, "Invalid global line in EIM potential file");
|
||||
}
|
||||
|
||||
division = values.next_double();
|
||||
rbig = values.next_double();
|
||||
rsmall = values.next_double();
|
||||
|
||||
found_global = true;
|
||||
} else if (type == "element:") {
|
||||
if (values.count() != 9) {
|
||||
error->one(FLERR, "Invalid element line in EIM potential file");
|
||||
}
|
||||
|
||||
std::string name = values.next_string();
|
||||
|
||||
ElementData data;
|
||||
data.ielement = values.next_int();
|
||||
data.mass = values.next_double();
|
||||
data.negativity = values.next_double();
|
||||
data.ra = values.next_double();
|
||||
data.ri = values.next_double();
|
||||
data.Ec = values.next_double();
|
||||
data.q0 = values.next_double();
|
||||
|
||||
if (elements.find(name) == elements.end()) {
|
||||
elements[name] = data;
|
||||
} else {
|
||||
error->one(FLERR, "Duplicate pair line in EIM potential file");
|
||||
}
|
||||
|
||||
} else if (type == "pair:") {
|
||||
if (values.count() != 17) {
|
||||
error->one(FLERR, "Invalid element line in EIM potential file");
|
||||
}
|
||||
|
||||
std::string elementA = values.next_string();
|
||||
std::string elementB = values.next_string();
|
||||
|
||||
PairData data;
|
||||
data.rcutphiA = values.next_double();
|
||||
data.rcutphiR = values.next_double();
|
||||
data.Eb = values.next_double();
|
||||
data.r0 = values.next_double();
|
||||
data.alpha = values.next_double();
|
||||
data.beta = values.next_double();
|
||||
data.rcutq = values.next_double();
|
||||
data.Asigma = values.next_double();
|
||||
data.rq = values.next_double();
|
||||
data.rcutsigma = values.next_double();
|
||||
data.Ac = values.next_double();
|
||||
data.zeta = values.next_double();
|
||||
data.rs = values.next_double();
|
||||
|
||||
// should be next_int, but since existing potential files have 1.0000e+00 format
|
||||
// we're doing this instead to keep compatibility
|
||||
data.tp = (int)values.next_double();
|
||||
|
||||
auto p = get_pair(elementA, elementB);
|
||||
|
||||
if (pairs.find(p) == pairs.end()) {
|
||||
pairs[p] = data;
|
||||
} else {
|
||||
error->one(FLERR, "Duplicate pair line in EIM potential file");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (!found_global) {
|
||||
error->one(FLERR, "Missing global line in EIM potential file");
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
ElementData & data = elements[name];
|
||||
setfl->ielement[i] = data.ielement;
|
||||
setfl->mass[i] = data.mass;
|
||||
setfl->negativity[i] = data.negativity;
|
||||
setfl->ra[i] = data.ra;
|
||||
setfl->ri[i] = data.ri;
|
||||
setfl->Ec[i] = data.Ec;
|
||||
setfl->q0[i] = data.q0;
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
|
||||
PairData & data = pairs[p];
|
||||
setfl->rcutphiA[ij] = data.rcutphiA;
|
||||
setfl->rcutphiR[ij] = data.rcutphiR;
|
||||
setfl->Eb[ij] = data.Eb;
|
||||
setfl->r0[ij] = data.r0;
|
||||
setfl->alpha[ij] = data.alpha;
|
||||
setfl->beta[ij] = data.beta;
|
||||
setfl->rcutq[ij] = data.rcutq;
|
||||
setfl->Asigma[ij] = data.Asigma;
|
||||
setfl->rq[ij] = data.rq;
|
||||
setfl->rcutsigma[ij] = data.rcutsigma;
|
||||
setfl->Ac[ij] = data.Ac;
|
||||
setfl->zeta[ij] = data.zeta;
|
||||
setfl->rs[ij] = data.rs;
|
||||
setfl->tp[ij] = data.tp;
|
||||
}
|
||||
|
||||
@ -21,6 +21,8 @@ PairStyle(eim,PairEIM)
|
||||
#define LMP_PAIR_EIM_H
|
||||
|
||||
#include "pair.h"
|
||||
#include <string>
|
||||
#include <map>
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
@ -40,16 +42,6 @@ class PairEIM : public Pair {
|
||||
void unpack_reverse_comm(int, int *, double *);
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
double **cutforcesq,cutmax;
|
||||
int nmax;
|
||||
double *rho,*fp;
|
||||
int rhofp;
|
||||
int *map; // which element each atom type maps to
|
||||
|
||||
int nelements; // # of elements to read from potential file
|
||||
char **elements; // element names
|
||||
|
||||
struct Setfl {
|
||||
double division,rbig,rsmall;
|
||||
int nr;
|
||||
@ -62,6 +54,17 @@ class PairEIM : public Pair {
|
||||
double ***Fij,***Gij,***phiij;
|
||||
double **cuts;
|
||||
};
|
||||
|
||||
protected:
|
||||
double **cutforcesq,cutmax;
|
||||
int nmax;
|
||||
double *rho,*fp;
|
||||
int rhofp;
|
||||
int *map; // which element each atom type maps to
|
||||
|
||||
int nelements; // # of elements to read from potential file
|
||||
char **elements; // element names
|
||||
|
||||
Setfl *setfl;
|
||||
|
||||
// potentials as array data
|
||||
@ -80,9 +83,6 @@ class PairEIM : public Pair {
|
||||
void allocate();
|
||||
void array2spline();
|
||||
void interpolate(int, double, double *, double **, double);
|
||||
int grabglobal(FILE *);
|
||||
int grabsingle(FILE *, int);
|
||||
int grabpair(FILE *, int, int);
|
||||
|
||||
double funccutoff(double, double, double);
|
||||
double funcphi(int, int, double);
|
||||
@ -94,6 +94,59 @@ class PairEIM : public Pair {
|
||||
void file2array();
|
||||
};
|
||||
|
||||
class EIMPotentialFileReader : protected Pointers {
|
||||
std::string filename;
|
||||
static const int MAXLINE = 1024;
|
||||
char line[MAXLINE];
|
||||
|
||||
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);
|
||||
|
||||
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
|
||||
double division;
|
||||
double rbig;
|
||||
double rsmall;
|
||||
|
||||
struct ElementData {
|
||||
int ielement;
|
||||
double mass;
|
||||
double negativity;
|
||||
double ra;
|
||||
double ri;
|
||||
double Ec;
|
||||
double q0;
|
||||
};
|
||||
|
||||
struct PairData {
|
||||
double rcutphiA;
|
||||
double rcutphiR;
|
||||
double Eb;
|
||||
double r0;
|
||||
double alpha;
|
||||
double beta;
|
||||
double rcutq;
|
||||
double Asigma;
|
||||
double rq;
|
||||
double rcutsigma;
|
||||
double Ac;
|
||||
double zeta;
|
||||
double rs;
|
||||
int tp;
|
||||
};
|
||||
|
||||
std::map<std::string, ElementData> elements;
|
||||
std::map<std::pair<std::string,std::string>, PairData> pairs;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
@ -2,3 +2,10 @@ add_executable(test_potential_file_reader test_potential_file_reader.cpp)
|
||||
target_link_libraries(test_potential_file_reader PRIVATE lammps GTest::GMock GTest::GTest)
|
||||
add_test(NAME PotentialFileReader COMMAND test_potential_file_reader WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
|
||||
set_tests_properties(PotentialFileReader PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
|
||||
|
||||
if (PKG_MANYBODY)
|
||||
add_executable(test_eim_potential_file_reader test_eim_potential_file_reader.cpp)
|
||||
target_link_libraries(test_eim_potential_file_reader PRIVATE lammps GTest::GMock GTest::GTest)
|
||||
add_test(NAME EIMPotentialFileReader COMMAND test_eim_potential_file_reader WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
|
||||
set_tests_properties(EIMPotentialFileReader PROPERTIES ENVIRONMENT "LAMMPS_POTENTIALS=${LAMMPS_POTENTIALS_DIR}")
|
||||
endif()
|
||||
|
||||
178
unittest/formats/test_eim_potential_file_reader.cpp
Normal file
178
unittest/formats/test_eim_potential_file_reader.cpp
Normal file
@ -0,0 +1,178 @@
|
||||
#include "gmock/gmock.h"
|
||||
#include "gtest/gtest.h"
|
||||
#include "lammps.h"
|
||||
#include "utils.h"
|
||||
#include "MANYBODY/pair_eim.h"
|
||||
|
||||
#include <mpi.h>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
class EIMPotentialFileReaderTest : public ::testing::Test {
|
||||
protected:
|
||||
LAMMPS * lmp;
|
||||
PairEIM::Setfl setfl;
|
||||
static const int nelements = 9;
|
||||
|
||||
void SetUp() override {
|
||||
const char *args[] = {"PotentialFileReaderTest", "-log", "none", "-echo", "screen", "-nocite" };
|
||||
char **argv = (char **)args;
|
||||
int argc = sizeof(args)/sizeof(char *);
|
||||
::testing::internal::CaptureStdout();
|
||||
lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
|
||||
::testing::internal::GetCapturedStdout();
|
||||
|
||||
int npair = nelements*(nelements+1)/2;
|
||||
setfl.ielement = new int[nelements];
|
||||
setfl.mass = new double[nelements];
|
||||
setfl.negativity = new double[nelements];
|
||||
setfl.ra = new double[nelements];
|
||||
setfl.ri = new double[nelements];
|
||||
setfl.Ec = new double[nelements];
|
||||
setfl.q0 = new double[nelements];
|
||||
setfl.rcutphiA = new double[npair];
|
||||
setfl.rcutphiR = new double[npair];
|
||||
setfl.Eb = new double[npair];
|
||||
setfl.r0 = new double[npair];
|
||||
setfl.alpha = new double[npair];
|
||||
setfl.beta = new double[npair];
|
||||
setfl.rcutq = new double[npair];
|
||||
setfl.Asigma = new double[npair];
|
||||
setfl.rq = new double[npair];
|
||||
setfl.rcutsigma = new double[npair];
|
||||
setfl.Ac = new double[npair];
|
||||
setfl.zeta = new double[npair];
|
||||
setfl.rs = new double[npair];
|
||||
setfl.tp = new int[npair];
|
||||
}
|
||||
|
||||
void TearDown() override {
|
||||
delete [] setfl.ielement;
|
||||
delete [] setfl.mass;
|
||||
delete [] setfl.negativity;
|
||||
delete [] setfl.ra;
|
||||
delete [] setfl.ri;
|
||||
delete [] setfl.Ec;
|
||||
delete [] setfl.q0;
|
||||
delete [] setfl.rcutphiA;
|
||||
delete [] setfl.rcutphiR;
|
||||
delete [] setfl.Eb;
|
||||
delete [] setfl.r0;
|
||||
delete [] setfl.alpha;
|
||||
delete [] setfl.beta;
|
||||
delete [] setfl.rcutq;
|
||||
delete [] setfl.Asigma;
|
||||
delete [] setfl.rq;
|
||||
delete [] setfl.rcutsigma;
|
||||
delete [] setfl.Ac;
|
||||
delete [] setfl.zeta;
|
||||
delete [] setfl.rs;
|
||||
delete [] setfl.tp;
|
||||
|
||||
::testing::internal::CaptureStdout();
|
||||
delete lmp;
|
||||
::testing::internal::GetCapturedStdout();
|
||||
}
|
||||
};
|
||||
|
||||
TEST_F(EIMPotentialFileReaderTest, global_line) {
|
||||
::testing::internal::CaptureStdout();
|
||||
EIMPotentialFileReader reader(lmp, "ffield.eim");
|
||||
::testing::internal::GetCapturedStdout();
|
||||
|
||||
reader.get_global(&setfl);
|
||||
ASSERT_DOUBLE_EQ(setfl.division, 2.0);
|
||||
ASSERT_DOUBLE_EQ(setfl.rbig, -1.645);
|
||||
ASSERT_DOUBLE_EQ(setfl.rsmall, 1.645);
|
||||
}
|
||||
|
||||
TEST_F(EIMPotentialFileReaderTest, element_line_sequential) {
|
||||
::testing::internal::CaptureStdout();
|
||||
EIMPotentialFileReader reader(lmp, "ffield.eim");
|
||||
::testing::internal::GetCapturedStdout();
|
||||
|
||||
reader.get_element(&setfl, 0, "Li");
|
||||
ASSERT_EQ(setfl.ielement[0], 3);
|
||||
ASSERT_DOUBLE_EQ(setfl.mass[0], 6.9410e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.negativity[0], 9.8000e-01);
|
||||
ASSERT_DOUBLE_EQ(setfl.ra[0], 1.1220e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.ri[0], 1.1220e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.Ec[0], -1.6500e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.q0[0], 0.0000e+00);
|
||||
|
||||
reader.get_element(&setfl, 1, "Na");
|
||||
ASSERT_EQ(setfl.ielement[1], 11);
|
||||
ASSERT_DOUBLE_EQ(setfl.mass[1], 2.2990e+01);
|
||||
ASSERT_DOUBLE_EQ(setfl.negativity[1], 9.3000e-01);
|
||||
ASSERT_DOUBLE_EQ(setfl.ra[1], 1.3690e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.ri[1], 1.3690e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.Ec[1], -1.1100e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.q0[1], 0.0000e+00);
|
||||
}
|
||||
|
||||
TEST_F(EIMPotentialFileReaderTest, element_line_random) {
|
||||
::testing::internal::CaptureStdout();
|
||||
EIMPotentialFileReader reader(lmp, "ffield.eim");
|
||||
::testing::internal::GetCapturedStdout();
|
||||
|
||||
reader.get_element(&setfl, 0, "Id");
|
||||
ASSERT_EQ(setfl.ielement[0], 53);
|
||||
ASSERT_DOUBLE_EQ(setfl.mass[0], 1.2690e+02);
|
||||
ASSERT_DOUBLE_EQ(setfl.negativity[0], 2.6600e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.ra[0], 1.8500e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.ri[0], 1.8500e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.Ec[0], -1.1100e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.q0[0], 0.0000e+00);
|
||||
}
|
||||
|
||||
TEST_F(EIMPotentialFileReaderTest, pair_line) {
|
||||
::testing::internal::CaptureStdout();
|
||||
EIMPotentialFileReader reader(lmp, "ffield.eim");
|
||||
::testing::internal::GetCapturedStdout();
|
||||
|
||||
reader.get_pair(&setfl, 0, "Li", "Li");
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutphiA[0], 6.0490e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutphiR[0], 6.0490e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.Eb[0], -2.5330e-01);
|
||||
ASSERT_DOUBLE_EQ(setfl.r0[0], 3.6176e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.alpha[0], 7.5536e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.beta[0], 3.5017e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutq[0], 0.0000e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.Asigma[0], 2.1778e-02);
|
||||
ASSERT_DOUBLE_EQ(setfl.rq[0], 2.0000e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutsigma[0], 7.0637e+00);
|
||||
ASSERT_DOUBLE_EQ(setfl.Ac[0], 3.3271e-01);
|
||||
ASSERT_DOUBLE_EQ(setfl.zeta[0], 6.0000e-01);
|
||||
ASSERT_DOUBLE_EQ(setfl.rs[0], 2.0000e+00);
|
||||
ASSERT_EQ(setfl.tp[0], 1);
|
||||
}
|
||||
|
||||
TEST_F(EIMPotentialFileReaderTest, pair_identical) {
|
||||
::testing::internal::CaptureStdout();
|
||||
EIMPotentialFileReader reader(lmp, "ffield.eim");
|
||||
::testing::internal::GetCapturedStdout();
|
||||
|
||||
reader.get_pair(&setfl, 0, "Li", "Na");
|
||||
reader.get_pair(&setfl, 1, "Na", "Li");
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutphiA[0], setfl.rcutphiA[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutphiR[0], setfl.rcutphiR[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.Eb[0], setfl.Eb[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.r0[0], setfl.r0[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.alpha[0], setfl.alpha[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.beta[0], setfl.beta[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutq[0], setfl.rcutq[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.Asigma[0], setfl.Asigma[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.rq[0], setfl.rq[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.rcutsigma[0], setfl.rcutsigma[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.Ac[0], setfl.Ac[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.zeta[0], setfl.zeta[1]);
|
||||
ASSERT_DOUBLE_EQ(setfl.rs[0], setfl.rs[1]);
|
||||
ASSERT_EQ(setfl.tp[0], setfl.tp[1]);
|
||||
}
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
MPI_Init(&argc, &argv);
|
||||
::testing::InitGoogleMock(&argc, argv);
|
||||
return RUN_ALL_TESTS();
|
||||
}
|
||||
Reference in New Issue
Block a user