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

update potential file parser in pair style extep

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-22 21:43:26 -04:00
parent 080662683d
commit 8ef9eae60a
4 changed files with 149 additions and 204 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
examples/PACKAGES/extep/BN.extep Symbolic link
View File

@ -0,0 +1 @@
../../../potentials/BN.extep

8
potentials/BN.extep
View File

@ -14,10 +14,10 @@ B B B 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498
N N N 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928 N N N 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928
B B N 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498959 2.5211820 2768.7363631 2.0 0.2 2.6857244 3376.3350735 B B N 3 1.0 0.0 26617.3000 141.2000 -0.1300 1.1422470 0.01498959 2.5211820 2768.7363631 2.0 0.2 2.6857244 3376.3350735
N N B 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928 N N B 3 1.0 0.0 23.5000 3.7500 -0.4000 0.6650000 0.01925100 2.6272721 2563.5603417 2.0 0.2 2.8293093 2978.9527928
B N B 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0 B N B 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2. 0.2 2.95 3330.0655849887
B N N 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0 B N N 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2. 0.2 2.95 3330.0655849887
N B B 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0 N B B 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2. 0.2 2.95 3330.0655849887
N B N 3 1.0 0.0d0 306.586555205d0 10.d0 -0.7218d0 0.6576543657d0 0.0027024851d0 2.69335d0 2595.6860833266d0 2.d0 0.2d0 2.95d0 3330.0655849887d0 N B N 3 1.0 0.0 306.586555205 10. -0.7218 0.6576543657 0.0027024851 2.69335 2595.6860833266 2. 0.2 2.95 3330.0655849887
# #
# 1.9925 Bicubic Splines Parameters # 1.9925 Bicubic Splines Parameters
# #

338
src/MANYBODY/pair_extep.cpp
View File

@ -28,6 +28,7 @@
#include "my_page.h" #include "my_page.h"
#include "neigh_list.h" #include "neigh_list.h"
#include "neighbor.h" #include "neighbor.h"
#include "potential_file_reader.h"
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
@ -511,234 +512,177 @@ double PairExTeP::init_one(int i, int j)
void PairExTeP::read_file(char *file) void PairExTeP::read_file(char *file)
{ {
int params_per_line = 17;
char **words = new char*[params_per_line+1];
memory->sfree(params); memory->sfree(params);
params = nullptr; params = nullptr;
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 = utils::open_potential(file,lmp,nullptr); PotentialFileReader reader(lmp, file, "ExTeP", unit_convert_flag);
if (fp == nullptr) char *line;
error->one(FLERR,"Cannot open ExTeP potential file {}: {}",file,utils::getsyserror());
}
// read each line out of file, skipping blank lines or leading '#' // transparently convert units for supported conversions
// store line of params if all 3 element tags are in element list
int n,nwords,ielement,jelement,kelement; int unit_convert = reader.get_unit_convert();
char line[MAXLINE],*ptr; double conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert);
int eof = 0;
while (true) { while ((line = reader.next_line(17))) {
if (comm->me == 0) { try {
ptr = fgets(line,MAXLINE,fp); ValueTokenizer values(line);
if (ptr == nullptr) {
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 std::string iname = values.next_string();
std::string jname = values.next_string();
std::string kname = values.next_string();
if ((ptr = strchr(line,'#'))) *ptr = '\0'; // ielement,jelement,kelement = 1st args
nwords = utils::count_words(line); // if all 3 args are in element list, then parse this line
if (nwords == 0) continue; // else skip to next entry in file
int ielement, jelement, kelement;
// concatenate additional lines until have params_per_line words 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;
while (nwords < params_per_line) { // load up parameter settings and error check their values
n = strlen(line);
if (comm->me == 0) { if (nparams == maxparam) {
ptr = fgets(&line[n],MAXLINE-n,fp); maxparam += DELTA;
if (ptr == nullptr) { params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
eof = 1;
fclose(fp); // make certain all addional allocated storage is initialized
} else n = strlen(line) + 1; // to avoid false positives when checking with valgrind
memset(params + nparams, 0, DELTA*sizeof(Param));
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].powerm = values.next_double();
params[nparams].gamma = values.next_double();
params[nparams].lam3 = values.next_double();
params[nparams].c = values.next_double();
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();
// currently only allow m exponent of 1 or 3
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());
} }
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break; if (params[nparams].c < 0.0 ||
MPI_Bcast(&n,1,MPI_INT,0,world); params[nparams].d < 0.0 ||
MPI_Bcast(line,n,MPI_CHAR,0,world); params[nparams].powern < 0.0 ||
if ((ptr = strchr(line,'#'))) *ptr = '\0'; params[nparams].beta < 0.0 ||
nwords = utils::count_words(line); 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].lam1 < 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)
error->one(FLERR,"Illegal ExTeP parameter");
nparams++;
if (nparams >= pow(nelements,3)) break;
} }
if (nwords != params_per_line) /* F_IJ (3) */
error->all(FLERR,"Insufficient spline parameters in potential file"); // initialize F_corr_data to all zeros
for (int iel=0; iel < nelements; iel++)
for (int jel=0; jel < nelements; jel++)
for (int in=0; in < 4; in++)
for (int jn=0; jn < 4; jn++)
for (int ivar=0; ivar < 3; ivar++)
F_corr_data[iel][jel][in][jn][ivar]=0;
// words = ptrs to all words in line // read the spline coefficients
nwords = 0; while ((line = reader.next_line(8))) {
words[nwords++] = strtok(line," \t\n\r\f"); try {
while ((words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue; ValueTokenizer values(line);
// ielement,jelement,kelement = 1st args std::string iname = values.next_string();
// if all 3 args are in element list, then parse this line std::string jname = values.next_string();
// else skip to next line std::string kname = values.next_string();
for (ielement = 0; ielement < nelements; ielement++) // ielement,jelement = 1st args
if (strcmp(words[0],elements[ielement]) == 0) break; // if all 2 args are in element list, then parse this line
if (ielement == nelements) continue; // else skip to next entry in file
for (jelement = 0; jelement < nelements; jelement++) int ielement, 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 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;
if (nparams == maxparam) { // skip line if it is a leftover from the previous section,
maxparam += DELTA; // which can be identified by having 3 elements (instead of 2)
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), // as first words.
"pair:params");
// make certain all addional allocated storage is initialized if (!utils::is_integer(kname))
// to avoid false positives when checking with valgrind continue;
memset(params + nparams, 0, DELTA*sizeof(Param)); int Ni = atoi(kname.c_str());
} int Nj = values.next_int();
double spline_val = values.next_double();
double spline_derx = values.next_double();
double spline_dery = values.next_double();
params[nparams].ielement = ielement; // Set value for all pairs of ielement,jelement (any kelement)
params[nparams].jelement = jelement; for (int iparam = 0; iparam < nparams; iparam++) {
params[nparams].kelement = kelement; if ( ielement == params[iparam].ielement && jelement == params[iparam].jelement) {
params[nparams].powerm = atof(words[3]); F_corr_data[ielement][jelement][Ni][Nj][0] = spline_val;
params[nparams].gamma = atof(words[4]); F_corr_data[ielement][jelement][Ni][Nj][1] = spline_derx;
params[nparams].lam3 = atof(words[5]); F_corr_data[ielement][jelement][Ni][Nj][2] = spline_dery;
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]);
// currently only allow m exponent of 1 or 3 F_corr_data[jelement][ielement][Nj][Ni][0] = spline_val;
F_corr_data[jelement][ielement][Nj][Ni][1] = spline_dery;
params[nparams].powermint = int(params[nparams].powerm); F_corr_data[jelement][ielement][Nj][Ni][2] = spline_derx;
}
if (params[nparams].c < 0.0 || params[nparams].d < 0.0 || }
params[nparams].powern < 0.0 || params[nparams].beta < 0.0 || } catch (TokenizerException &e) {
params[nparams].lam2 < 0.0 || params[nparams].bigb < 0.0 || error->one(FLERR, e.what());
params[nparams].bigr < 0.0 ||params[nparams].bigd < 0.0 ||
params[nparams].bigd > params[nparams].bigr ||
params[nparams].lam1 < 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)
error->all(FLERR,"Illegal ExTeP parameter");
nparams++;
if (nparams >= pow(nelements,3)) break;
}
// deallocate words array
delete[] words;
/* F_IJ (3) */
// read the spline coefficients
params_per_line = 8;
// reallocate with new size
words = new char*[params_per_line+1];
// initialize F_corr_data to all zeros
for (int iel=0;iel<nelements;iel++)
for (int jel=0;jel<nelements;jel++)
for (int in=0;in<4;in++)
for (int jn=0;jn<4;jn++)
for (int ivar=0;ivar<3;ivar++)
F_corr_data[iel][jel][in][jn][ivar]=0;
// loop until EOF
while (true) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
//fputs(line,stdout);
if (ptr == nullptr) {
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
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = utils::count_words(line);
if (nwords == 0) continue;
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((nwords < params_per_line)
&& (words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue;
// skip line if it is a leftover from the previous section,
// which can be identified by having 3 elements (instead of 2)
// as first words.
if (isupper(words[0][0]) && isupper(words[1][0]) && isupper(words[2][0]))
continue;
// need to have two elements followed by a number in each line
if (!(isupper(words[0][0]) && isupper(words[1][0])
&& !isupper(words[2][0])))
error->all(FLERR,"Incorrect format in ExTeP potential file");
// ielement,jelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next line
// these lines set ielement and jelement to the
// integers matching the strings from the input
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;
int Ni = atoi(words[2]);
int Nj = atoi(words[3]);
double spline_val = atof(words[4]);
double spline_derx = atof(words[5]);
double spline_dery = atof(words[6]);
// Set value for all pairs of ielement,jelement (any kelement)
for (int iparam = 0; iparam < nparams; iparam++) {
if ( ielement == params[iparam].ielement
&& jelement == params[iparam].jelement) {
F_corr_data[ielement][jelement][Ni][Nj][0] = spline_val;
F_corr_data[ielement][jelement][Ni][Nj][1] = spline_derx;
F_corr_data[ielement][jelement][Ni][Nj][2] = spline_dery;
F_corr_data[jelement][ielement][Nj][Ni][0] = spline_val;
F_corr_data[jelement][ielement][Nj][Ni][1] = spline_dery;
F_corr_data[jelement][ielement][Nj][Ni][2] = spline_derx;
} }
} }
} }
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
delete[] words; if (comm->me != 0) {
/* END F_IJ (3) */ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
MPI_Bcast(&F_corr_data[0][0][0][0][0], MAXTYPES*MAXTYPES*NSPLINE*NSPLINE*3, MPI_DOUBLE, 0, world);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

6
src/MANYBODY/pair_tersoff.cpp
View File

@ -428,6 +428,7 @@ void PairTersoff::read_file(char *file)
int unit_convert = reader.get_unit_convert(); int unit_convert = reader.get_unit_convert();
double conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert); double conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert);
while ((line = reader.next_line(NPARAMS_PER_LINE))) { while ((line = reader.next_line(NPARAMS_PER_LINE))) {
try { try {
ValueTokenizer values(line); ValueTokenizer values(line);
@ -455,8 +456,7 @@ void PairTersoff::read_file(char *file)
if (nparams == maxparam) { if (nparams == maxparam) {
maxparam += DELTA; maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
"pair:params");
// make certain all addional allocated storage is initialized // make certain all addional allocated storage is initialized
// to avoid false positives when checking with valgrind // to avoid false positives when checking with valgrind
@ -505,7 +505,7 @@ void PairTersoff::read_file(char *file)
params[nparams].biga < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].powerm - params[nparams].powermint != 0.0 || params[nparams].powerm - params[nparams].powermint != 0.0 ||
(params[nparams].powermint != 3 && (params[nparams].powermint != 3 &&
params[nparams].powermint != 1) || params[nparams].powermint != 1) ||
params[nparams].gamma < 0.0) params[nparams].gamma < 0.0)
error->one(FLERR,"Illegal Tersoff parameter"); error->one(FLERR,"Illegal Tersoff parameter");