From 14fc42833f02832c611871b9fcc9ed8a7a101bd3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 10 Dec 2021 08:45:01 -0500 Subject: [PATCH] modernize potential file reader for local/density --- src/MANYBODY/pair_local_density.cpp | 132 +++++++++++++++------------- 1 file changed, 71 insertions(+), 61 deletions(-) diff --git a/src/MANYBODY/pair_local_density.cpp b/src/MANYBODY/pair_local_density.cpp index 26b24d8445..3fe12bb71d 100644 --- a/src/MANYBODY/pair_local_density.cpp +++ b/src/MANYBODY/pair_local_density.cpp @@ -28,6 +28,7 @@ #include "memory.h" #include "neigh_list.h" #include "neighbor.h" +#include "potential_file_reader.h" #include @@ -656,40 +657,37 @@ void PairLocalDensity::interpolate_cbspl(int n, double delta, void PairLocalDensity::parse_file(char *filename) { - int k, n; - int me = comm->me; - FILE *fptr; - char line[MAXLINE]; - double ratio, lc2, uc2, denom; + // parse potential file header + if (comm->me == 0) { + PotentialFileReader reader(lmp, filename, "local/density"); - if (me == 0) { - fptr = fopen(filename, "r"); - if (fptr == nullptr) - error->one(FLERR,"Cannot open Local Density potential file {}: {}",filename,utils::getsyserror()); + try { + + // ignore first 2 comment lines + reader.skip_line(); + reader.skip_line(); + + // extract number of potentials and number of (frho, rho) points + ValueTokenizer values = reader.next_values(2); + nLD = values.next_int(); + nrho = values.next_int(); + + const int numld = atom->ntypes*atom->ntypes; + if (nLD != numld) + error->warning(FLERR, "Expected {} local density potentials but got {}",numld, nLD); + + } catch (TokenizerException &e) { + error->one(FLERR, e.what()); + } } - double *ftmp; // tmp var to extract the complete 2D frho array from file - - // broadcast number of LD potentials and number of (rho,frho) pairs - if (me == 0) { - - // first 2 comment lines ignored - utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error); - utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error); - - // extract number of potentials and number of (frho, rho) points - utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error); - sscanf(line, "%d %d", &nLD, &nrho); - utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error); - } + // broadcast number of LD potentials and number of (rho,frho) pairs and allocate storage MPI_Bcast(&nLD,1,MPI_INT,0,world); MPI_Bcast(&nrho,1,MPI_INT,0,world); comm_forward = comm_reverse = nLD; - if ((me == 0) && (nLD != atom->ntypes*atom->ntypes)) - error->warning(FLERR, "Expected {} local density potentials but got {}", - atom->ntypes*atom->ntypes, nLD); + double *ftmp; // tmp var to extract the complete 2D frho array from file // setting up all arrays to be read from files and broadcasted memory->create(uppercut, nLD, "pairLD:uppercut"); @@ -708,54 +706,65 @@ void PairLocalDensity::parse_file(char *filename) { // setting up central and neighbor atom filters memory->create(a, nLD, atom->ntypes+1 , "pairLD:a"); memory->create(b, nLD, atom->ntypes+1, "pairLD:b"); - if (me == 0) { - for (n = 1; n <= atom->ntypes; n++) { - for (k = 0; k < nLD; k++) { - a[k][n] = 0; - b[k][n] = 0; - } + for (int k = 0; k < nLD; k++) { + for (int n = 1; n <= atom->ntypes; n++) { + a[k][n] = 0; + b[k][n] = 0; } } - // read file block by block + // parse potential file body + if (comm->me == 0) { + PotentialFileReader reader(lmp, filename, "local/density"); - if (me == 0) { - for (k = 0; k < nLD; k++) { + try { + double ratio, lc2, uc2, denom; + ValueTokenizer values(""); + + // ignore first 4 lines already processed + + reader.skip_line(); + reader.skip_line(); + reader.skip_line(); + reader.skip_line(); + + for (int k = 0; k < nLD; k++) { // parse upper and lower cut values - if (fgets(line,MAXLINE,fptr)==nullptr) break; - sscanf(line, "%lf %lf", &lowercut[k], &uppercut[k]); + values = reader.next_values(2); + lowercut[k] = values.next_double(); + uppercut[k] = values.next_double(); - // parse and broadcast central atom filter - utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error); - char *tmp = strtok(line, " /t/n/r/f"); - while (tmp != nullptr) { - a[k][atoi(tmp)] = 1; - tmp = strtok(nullptr, " /t/n/r/f"); + // parse central atom filter + values = ValueTokenizer(reader.next_line()); + while (values.has_next()) { + int atype = values.next_int(); + if ((atype < 1) || (atype > atom->ntypes)) + throw TokenizerException("Invalid atom type filter value",std::to_string(atype)); + a[k][atype] = 1; } // parse neighbor atom filter - utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error); - tmp = strtok(line, " /t/n/r/f"); - while (tmp != nullptr) { - b[k][atoi(tmp)] = 1; - tmp = strtok(nullptr, " /t/n/r/f"); + values = ValueTokenizer(reader.next_line()); + while (values.has_next()) { + int btype = values.next_int(); + if ((btype < 1) || (btype > atom->ntypes)) + throw TokenizerException("Invalid atom type filter value",std::to_string(btype)); + b[k][btype] = 1; } // parse min, max and delta rho values - utils::sfgets(FLERR,line, MAXLINE, fptr,filename,error); - sscanf(line, "%lf %lf %lf", &rho_min[k], &rho_max[k], &delta_rho[k]); + values = reader.next_values(3); + rho_min[k] = values.next_double(); + rho_max[k] = values.next_double(); // recompute delta_rho from scratch for precision delta_rho[k] = (rho_max[k] - rho_min[k]) / (nrho - 1); // parse tabulated frho values from each line into temporary array - for (n = 0; n < nrho; n++) { - utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error); - sscanf(line, "%lf", &ftmp[k*nrho + n]); - } + reader.next_dvector(ftmp+k*nrho, nrho); // ignore blank line at the end of every block - utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error); + reader.skip_line(); // set coefficients for local density indicator function uc2 = uppercut[k] * uppercut[k]; @@ -770,6 +779,10 @@ void PairLocalDensity::parse_file(char *filename) { c4[k] = -(3.0 + 3.0*ratio) / (uc2*uc2 * denom); c6[k] = 2.0 / (uc2*uc2*uc2 * denom); } + + } catch (TokenizerException &e) { + error->one(FLERR, e.what()); + } } // Broadcast all parsed arrays @@ -781,7 +794,7 @@ void PairLocalDensity::parse_file(char *filename) { MPI_Bcast(&c2[0], nLD, MPI_DOUBLE, 0, world); MPI_Bcast(&c4[0], nLD, MPI_DOUBLE, 0, world); MPI_Bcast(&c6[0], nLD, MPI_DOUBLE, 0, world); - for (k = 0; k < nLD; k++) { + for (int k = 0; k < nLD; k++) { MPI_Bcast(&a[k][1], atom->ntypes, MPI_INT, 0, world); MPI_Bcast(&b[k][1], atom->ntypes, MPI_INT, 0, world); } @@ -790,14 +803,12 @@ void PairLocalDensity::parse_file(char *filename) { MPI_Bcast(&delta_rho[0], nLD, MPI_DOUBLE, 0, world); MPI_Bcast(&ftmp[0], nLD*nrho, MPI_DOUBLE, 0, world); - if (me == 0) fclose(fptr); - // set up rho and frho arrays memory->create(rho, nLD, nrho, "pairLD:rho"); memory->create(frho, nLD, nrho, "pairLD:frho"); - for (k = 0; k < nLD; k++) { - for (n = 0; n < nrho; n++) { + for (int k = 0; k < nLD; k++) { + for (int n = 0; n < nrho; n++) { rho[k][n] = rho_min[k] + n*delta_rho[k]; frho[k][n] = ftmp[k*nrho + n]; } @@ -885,4 +896,3 @@ double PairLocalDensity::memory_usage() bytes += (double)2 * (nmax*nLD) * sizeof(double); return bytes; } -