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

refactor file reader code

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-16 15:29:43 -04:00
parent 9e559de601
commit 436d200a9a
2 changed files with 74 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files

180
src/INTERLAYER/pair_kolmogorov_crespi_full.cpp
View File

@ -33,6 +33,8 @@
#include "neigh_list.h" #include "neigh_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "neighbor.h" #include "neighbor.h"
#include "potential_file_reader.h"
#include "tokenizer.h"
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
@ -183,138 +185,103 @@ double PairKolmogorovCrespiFull::init_one(int i, int j)
void PairKolmogorovCrespiFull::read_file(char *filename) void PairKolmogorovCrespiFull::read_file(char *filename)
{ {
int params_per_line = 12;
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(filename,lmp,nullptr); PotentialFileReader reader(lmp, filename, "kolmogorov/crespi/full", unit_convert_flag);
if (fp == nullptr) char *line;
error->one(FLERR,"Cannot open KC potential file {}: {}",filename,utils::getsyserror());
}
// 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
int i,j,n,m,nwords,ielement,jelement; try {
char line[MAXLINE],*ptr; ValueTokenizer values(line);
int eof = 0;
while (1) { std::string iname = values.next_string();
if (comm->me == 0) { std::string jname = values.next_string();
ptr = fgets(line,MAXLINE,fp);
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 // ielement,jelement = 1st args
// if both args are in element list, then parse this line
// else skip to next entry in file
int ielement, jelement;
if ((ptr = strchr(line,'#'))) *ptr = '\0'; for (ielement = 0; ielement < nelements; ielement++)
nwords = utils::count_words(line); if (iname == elements[ielement]) break;
if (nwords == 0) continue; if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (jname == elements[jelement]) break;
if (jelement == nelements) continue;
// concatenate additional lines until have params_per_line words // expand storage, if needed
while (nwords < params_per_line) { if (nparams == maxparam) {
n = strlen(line); maxparam += DELTA;
if (comm->me == 0) { params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == nullptr) { // make certain all addional allocated storage is initialized
eof = 1; // to avoid false positives when checking with valgrind
fclose(fp);
} else n = strlen(line) + 1; memset(params + nparams, 0, DELTA*sizeof(Param));
}
// load up parameter settings and error check their values
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].z0 = values.next_double();
params[nparams].C0 = values.next_double();
params[nparams].C2 = values.next_double();
params[nparams].C4 = values.next_double();
params[nparams].C = values.next_double();
params[nparams].delta = values.next_double();
params[nparams].lambda = values.next_double();
params[nparams].A = values.next_double();
// S provides a convenient scaling of all energies
params[nparams].S = values.next_double();
params[nparams].rcut = values.next_double();
} catch (TokenizerException &e) {
error->one(FLERR, e.what());
} }
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break; // energies in meV further scaled by S
MPI_Bcast(&n,1,MPI_INT,0,world); const double meV = 1.0e-3*params[nparams].S;
MPI_Bcast(line,n,MPI_CHAR,0,world); params[nparams].C *= meV;
if ((ptr = strchr(line,'#'))) *ptr = '\0'; params[nparams].A *= meV;
nwords = utils::count_words(line); params[nparams].C0 *= meV;
params[nparams].C2 *= meV;
params[nparams].C4 *= meV;
// precompute some quantities
params[nparams].delta2inv = pow(params[nparams].delta,-2);
params[nparams].z06 = pow(params[nparams].z0,6);
nparams++;
} }
if (nwords != params_per_line) MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
error->all(FLERR,"Insufficient format in KC potential file"); MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
// words = ptrs to all words in line if (comm->me != 0) {
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue;
// ielement,jelement = 1st args
// if these 2 args are in element list, then parse this line
// else skip to next line (continue)
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;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
// make certain all addional allocated storage is initialized
// to avoid false positives when checking with valgrind
memset(params + nparams, 0, DELTA*sizeof(Param));
} }
params[nparams].ielement = ielement; MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
params[nparams].jelement = jelement;
params[nparams].z0 = atof(words[2]);
params[nparams].C0 = atof(words[3]);
params[nparams].C2 = atof(words[4]);
params[nparams].C4 = atof(words[5]);
params[nparams].C = atof(words[6]);
params[nparams].delta = atof(words[7]);
params[nparams].lambda = atof(words[8]);
params[nparams].A = atof(words[9]);
// S provides a convenient scaling of all energies
params[nparams].S = atof(words[10]);
params[nparams].rcut = atof(words[11]);
// energies in meV further scaled by S
double meV = 1.0e-3*params[nparams].S;
params[nparams].C *= meV;
params[nparams].A *= meV;
params[nparams].C0 *= meV;
params[nparams].C2 *= meV;
params[nparams].C4 *= meV;
// precompute some quantities
params[nparams].delta2inv = pow(params[nparams].delta,-2);
params[nparams].z06 = pow(params[nparams].z0,6);
nparams++;
//if(nparams >= pow(atom->ntypes,3)) break;
} }
memory->destroy(elem2param); memory->destroy(elem2param);
memory->destroy(cutKCsq); memory->destroy(cutKCsq);
memory->create(elem2param,nelements,nelements,"pair:elem2param"); memory->create(elem2param,nelements,nelements,"pair:elem2param");
memory->create(cutKCsq,nelements,nelements,"pair:cutKCsq"); memory->create(cutKCsq,nelements,nelements,"pair:cutKCsq");
for (i = 0; i < nelements; i++) { for (int i = 0; i < nelements; i++) {
for (j = 0; j < nelements; j++) { for (int j = 0; j < nelements; j++) {
n = -1; int n = -1;
for (m = 0; m < nparams; m++) { for (int m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement) { if (i == params[m].ielement && j == params[m].jelement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry"); if (n >= 0) error->all(FLERR,"KC Potential file has duplicate entry");
n = m; n = m;
} }
} }
@ -323,7 +290,6 @@ void PairKolmogorovCrespiFull::read_file(char *filename)
cutKCsq[i][j] = params[n].rcut*params[n].rcut; cutKCsq[i][j] = params[n].rcut*params[n].rcut;
} }
} }
delete [] words;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

1
src/INTERLAYER/pair_kolmogorov_crespi_full.h
View File

@ -39,6 +39,7 @@ class PairKolmogorovCrespiFull : public Pair {
void calc_FRep(int, int); void calc_FRep(int, int);
void calc_FvdW(int, int); void calc_FvdW(int, int);
double single(int, int, int, int, double, double, double, double &); double single(int, int, int, int, double, double, double, double &);
static constexpr int NPARAMS_PER_LINE = 12;
protected: protected:
int me; int me;