From 38bcc493b011ca9bc4e1b47efc7fd57870dc73c8 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 10:28:40 -0500 Subject: [PATCH 1/8] simplify and modernize CMAP file parser --- src/MOLECULE/fix_cmap.cpp | 148 +++----------------------------------- 1 file changed, 11 insertions(+), 137 deletions(-) diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index a763c5d14c..c6e356e829 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -37,6 +37,7 @@ #include "force.h" #include "math_const.h" #include "memory.h" +#include "potential_file_reader.h" #include "respa.h" #include "tokenizer.h" #include "update.h" @@ -634,150 +635,23 @@ double FixCMAP::compute_scalar() void FixCMAP::read_grid_map(char *cmapfile) { - char linebuf[MAXLINE]; - char *chunk,*line; - int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter; - - FILE *fp = nullptr; if (comm->me == 0) { - fp = utils::open_potential(cmapfile,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open fix cmap file {}: {}", - cmapfile, utils::getsyserror()); + try { + memset(&cmapgrid[0][0][0], 6*CMAPDIM*CMAPDIM, sizeof(double)); + PotentialFileReader reader(lmp, cmapfile, "cmap grid"); - } + // there are six maps in this order. + // alanine, alanine-proline, proline, proline-proline, glycine, glycine-proline. + // read as one big blob of numbers while ignoring comments - for (int ix1 = 0; ix1 < 6; ix1++) - for (int ix2 = 0; ix2 < CMAPDIM; ix2++) - for (int ix3 = 0; ix3 < CMAPDIM; ix3++) - cmapgrid[ix1][ix2][ix3] = 0.0; + reader.next_dvector(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM); - counter = 0; - i1 = i2 = i3 = i4 = i5 = i6 = 0; - j1 = j2 = j3 = j4 = j5 = j6 = 0; - - int done = 0; - - while (!done) { - // only read on rank 0 and broadcast to all other ranks - if (comm->me == 0) - done = (fgets(linebuf,MAXLINE,fp) == nullptr); - - MPI_Bcast(&done,1,MPI_INT,0,world); - if (done) continue; - - MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world); - - // remove leading whitespace - line = linebuf; - while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line; - - // skip if empty line or comment - if (!line || *line =='\n' || *line == '\0' || *line == '#') continue; - - // read in the cmap grid point values - // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors - // will occur if content of the file "cmap.data" is altered - // - // Reading order of the maps: - // 1. Alanine map - // 2. Alanine before proline map - // 3. Proline map - // 4. Two adjacent prolines map - // 5. Glycine map - // 6. Glycine before proline map - - chunk = strtok(line, " \r\n"); - while (chunk != nullptr) { - - // alanine map - - if (counter < CMAPDIM*CMAPDIM) { - cmapgrid[0][i1][j1] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j1++; - if (j1 == CMAPDIM) { - j1 = 0; - i1++; - } - counter++; - } - - // alanine-proline map - - else if (counter >= CMAPDIM*CMAPDIM && - counter < 2*CMAPDIM*CMAPDIM) { - cmapgrid[1][i2][j2]= atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j2++; - if (j2 == CMAPDIM) { - j2 = 0; - i2++; - } - counter++; - } - - // proline map - - else if (counter >= 2*CMAPDIM*CMAPDIM && - counter < 3*CMAPDIM*CMAPDIM) { - cmapgrid[2][i3][j3] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j3++; - if (j3 == CMAPDIM) { - j3 = 0; - i3++; - } - counter++; - } - - // 2 adjacent prolines map - - else if (counter >= 3*CMAPDIM*CMAPDIM && - counter < 4*CMAPDIM*CMAPDIM) { - cmapgrid[3][i4][j4] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j4++; - if (j4 == CMAPDIM) { - j4 = 0; - i4++; - } - counter++; - } - - // glycine map - - else if (counter >= 4*CMAPDIM*CMAPDIM && - counter < 5*CMAPDIM*CMAPDIM) { - cmapgrid[4][i5][j5] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j5++; - if (j5 == CMAPDIM) { - j5 = 0; - i5++; - } - counter++; - } - - // glycine-proline map - - else if (counter >= 5*CMAPDIM*CMAPDIM && - counter < 6*CMAPDIM*CMAPDIM) { - cmapgrid[5][i6][j6] = atof(chunk); - chunk = strtok(nullptr, " \r\n"); - j6++; - if (j6 == CMAPDIM) { - j6 = 0; - i6++; - } - counter++; - } - - else break; + } catch (std::exception &e) { + error->one(FLERR,"Error reading CMAP potential file: {}", e.what()); } } - if (comm->me == 0) fclose(fp); + MPI_Bcast(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ From ef13455d6b3ad60cf494b79a588a74a8fb8a4f67 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 16:48:07 -0500 Subject: [PATCH 2/8] add some metadata tags --- examples/PACKAGES/electron_stopping/Si.Si.elstop | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/PACKAGES/electron_stopping/Si.Si.elstop b/examples/PACKAGES/electron_stopping/Si.Si.elstop index 51618149ff..1a6952e571 100644 --- a/examples/PACKAGES/electron_stopping/Si.Si.elstop +++ b/examples/PACKAGES/electron_stopping/Si.Si.elstop @@ -1,8 +1,8 @@ -# Electronic stopping for Si in Si -# Uses metal units +# Electronic stopping for Si in Si UNITS: metal DATE: 2019-02-19 CONTRIBUTOR: Risto Toijala risto.toijala@helsinki.fi # Kinetic energy in eV, stopping power in eV/A # For other atom types, add columns. +# atom type 1 # energy Si in Si 3918.2 6.541 15672.9 13.091 From b890b564caa7571536c750ab16e148018e22fb03 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 16:48:26 -0500 Subject: [PATCH 3/8] modernize electron stopping file reader --- src/EXTRA-FIX/fix_electron_stopping.cpp | 57 +++++++++++-------------- 1 file changed, 26 insertions(+), 31 deletions(-) diff --git a/src/EXTRA-FIX/fix_electron_stopping.cpp b/src/EXTRA-FIX/fix_electron_stopping.cpp index bd54fea97d..49fcdd41bf 100644 --- a/src/EXTRA-FIX/fix_electron_stopping.cpp +++ b/src/EXTRA-FIX/fix_electron_stopping.cpp @@ -29,6 +29,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include "region.h" #include "update.h" @@ -236,49 +237,43 @@ double FixElectronStopping::compute_scalar() return SeLoss_all; } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + read electron stopping parameters. only called from MPI rank 0. + format: energy then one column per atom type + read as many lines as available. + energies must be sorted in ascending order. + ---------------------------------------------------------------------- */ void FixElectronStopping::read_table(const char *file) { - char line[MAXLINE]; - - FILE *fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open stopping range table {}: {}", file, utils::getsyserror()); - const int ncol = atom->ntypes + 1; + int nlines = 0; + PotentialFileReader reader(lmp, file, "electron stopping data table"); - int l = 0; - while (true) { - if (fgets(line, MAXLINE, fp) == nullptr) break; // end of file - if (line[0] == '#') continue; // comment + try { + char *line; + double oldvalue = 0.0; - char *pch = strtok(line, " \t\n\r"); - if (pch == nullptr) continue; // blank line + while ((line = reader.next_line())) { + if (nlines >= maxlines) grow_table(); + ValueTokenizer values(line); + elstop_ranges[0][nlines] = values.next_double(); + if (elstop_ranges[0][nlines] <= oldvalue) + throw TokenizerException("energy values must be positive and in ascending order",line); - if (l >= maxlines) grow_table(); + oldvalue = elstop_ranges[0][nlines]; + for (int i = 1; i < ncol; ++i) + elstop_ranges[i][nlines] = values.next_double(); - int i = 0; - for ( ; i < ncol && pch != nullptr; i++) { - elstop_ranges[i][l] = utils::numeric(FLERR, pch,false,lmp); - pch = strtok(nullptr, " \t\n\r"); + ++nlines; } - - if (i != ncol || pch != nullptr) // too short or too long - error->one(FLERR, "fix electron/stopping: Invalid table line"); - - if (l >= 1 && elstop_ranges[0][l] <= elstop_ranges[0][l-1]) - error->one(FLERR, - "fix electron/stopping: Energies must be in ascending order"); - - l++; + } catch (std::exception &e) { + error->one(FLERR, "Problem parsing electron stopping data: {}", e.what()); } - table_entries = l; - - if (table_entries == 0) + if (nlines == 0) error->one(FLERR, "Did not find any data in electron/stopping table file"); - fclose(fp); + table_entries = nlines; } /* ---------------------------------------------------------------------- */ From e22f62f6db5ce02e8571150eab62437323e907d4 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 19:55:30 -0500 Subject: [PATCH 4/8] add some metadata for fix gle matrix files --- examples/PACKAGES/gle/qt-300k.A | 2 ++ examples/PACKAGES/gle/qt-300k.C | 2 ++ examples/PACKAGES/gle/smart.A | 2 ++ 3 files changed, 6 insertions(+) diff --git a/examples/PACKAGES/gle/qt-300k.A b/examples/PACKAGES/gle/qt-300k.A index 073100ad9d..31ada37b05 100644 --- a/examples/PACKAGES/gle/qt-300k.A +++ b/examples/PACKAGES/gle/qt-300k.A @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 3.904138445158e-4 4.681059722010e-2 3.079778738058e-2 4.428079381336e-2 5.057825203477e-2 2.591193419597e-2 1.513907125942e-2 -4.789343294190e-2 2.037551040972e-2 6.597801861779e-2 -8.528273506602e-3 -2.230839572773e-3 6.593086307870e-3 -6.653653981891e-2 -2.905096373618e-2 -6.597801861779e-2 2.086885297843e-2 1.145471984072e-2 3.111465343867e-2 1.101562087523e-2 -3.264959166486e-2 diff --git a/examples/PACKAGES/gle/qt-300k.C b/examples/PACKAGES/gle/qt-300k.C index 50d68df251..d6239ca8dc 100644 --- a/examples/PACKAGES/gle/qt-300k.C +++ b/examples/PACKAGES/gle/qt-300k.C @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 2.999985914100e+2 5.937593337000e+0 2.144376751500e+2 5.883055908000e+1 -1.119803766000e+2 -6.793381764000e+1 1.379789732400e+1 5.937593337000e+0 3.781851303000e+2 -5.794518522000e+1 -2.178772681500e+2 -1.649310659100e+2 -6.557113452000e+1 3.833830743000e+1 2.144376751500e+2 -5.794518522000e+1 7.325759985000e+2 2.188507713000e+2 -3.704586531000e+2 -3.385193865000e+1 -4.827706758000e+0 diff --git a/examples/PACKAGES/gle/smart.A b/examples/PACKAGES/gle/smart.A index a63bd881c8..8abbdea76a 100644 --- a/examples/PACKAGES/gle/smart.A +++ b/examples/PACKAGES/gle/smart.A @@ -1,3 +1,5 @@ +# DATE: 2014-11-24 UNITS: real CONTRIBUTOR: Michele Ceriotti michele.ceriotti@gmail.com + 4.247358737218e-3 3.231404593065e-2 -9.715629522215e-3 8.199488441225e-3 7.288427565896e-3 -3.634229949055e-3 -3.294395982200e-3 4.887738278128e-2 5.978602802893e-1 -2.079222967202e-1 1.088740438247e-1 4.666954324786e-2 -1.334147134493e-2 -3.914996645811e-2 4.015863091190e-2 2.079222967202e-1 1.645970119444e-1 8.351952991453e-2 5.114740085468e-2 1.217862237677e-2 -4.506711205974e-2 From c4c3b545b2f4c3507dbeb2696c5f8e062a1af037 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 18 Jan 2022 19:55:50 -0500 Subject: [PATCH 5/8] modernize fix gle matrix file reading --- src/EXTRA-FIX/fix_gle.cpp | 104 ++++++++++---------------------------- 1 file changed, 26 insertions(+), 78 deletions(-) diff --git a/src/EXTRA-FIX/fix_gle.cpp b/src/EXTRA-FIX/fix_gle.cpp index de225c6980..d7061e199c 100644 --- a/src/EXTRA-FIX/fix_gle.cpp +++ b/src/EXTRA-FIX/fix_gle.cpp @@ -14,7 +14,7 @@ /* ---------------------------------------------------------------------- Contributing authors: Michele Ceriotti (EPFL), Joe Morrone (Stony Brook), - Axel Kohylmeyer (Temple U) + Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ #include "fix_gle.h" @@ -24,6 +24,7 @@ #include "error.h" #include "force.h" #include "memory.h" +#include "potential_file_reader.h" #include "random_mars.h" #include "respa.h" #include "update.h" @@ -221,46 +222,17 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : int seed = utils::inumeric(FLERR,arg[6],false,lmp); // LOADING A matrix - FILE *fgle = nullptr; char *fname = arg[7]; + memset(A, ns1sq, sizeof(double)); if (comm->me == 0) { - fgle = utils::open_potential(fname,lmp,nullptr); - if (fgle == nullptr) - error->one(FLERR,"Cannot open A-matrix file {}: {}",fname, utils::getsyserror()); - utils::logmesg(lmp,"Reading A-matrix from {}\n", fname); - } - - // read each line of the file, skipping blank lines or leading '#' - - char line[MAXLINE],*ptr; - int n,nwords,ndone=0,eof=0; - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fgle); - if (ptr == nullptr) { - eof = 1; - fclose(fgle); - } else n = strlen(line) + 1; + PotentialFileReader reader(lmp,fname,"fix gle A matrix"); + try { + reader.next_dvector(A, ns1sq); + } catch (std::exception &e) { + error->one(FLERR,"Error reading A-matrix: {}", e.what()); } - 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; - - ptr = strtok(line," \t\n\r\f"); - do { - A[ndone] = atof(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - ndone++; - } while ((ptr != nullptr) && (ndone < ns1sq)); } + MPI_Bcast(A,ns1sq,MPI_DOUBLE,0,world); fnoneq=0; gle_every=1; gle_step=0; for (int iarg=8; iargme == 0) { - fgle = utils::open_potential(fname,lmp,nullptr); - if (fgle == nullptr) - error->one(FLERR,"Cannot open C-matrix file {}: {}",fname, utils::getsyserror()); - utils::logmesg(lmp,"Reading C-matrix from {}\n", fname); - } - - // read each line of the file, skipping blank lines or leading '#' - ndone = eof = 0; - const double cfac = force->boltz / force->mvv2e; - - while (true) { - if (comm->me == 0) { - ptr = fgets(line,MAXLINE,fgle); - if (ptr == nullptr) { - eof = 1; - fclose(fgle); - } else n = strlen(line) + 1; + PotentialFileReader reader(lmp,fname,"fix gle C matrix"); + try { + reader.next_dvector(C, ns1sq); + } catch (std::exception &e) { + error->one(FLERR,"Error reading C-matrix: {}", e.what()); } - 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; - - ptr = strtok(line," \t\n\r\f"); - do { - C[ndone] = cfac*atof(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - ndone++; - } while ((ptr != nullptr) && (ndone < ns1sq)); } + MPI_Bcast(C,ns1sq,MPI_DOUBLE,0,world); + const double cfac = force->boltz / force->mvv2e; + for (int i=0; i < ns1sq; ++i) C[i] *= cfac; } #ifdef GLE_DEBUG @@ -366,12 +314,12 @@ FixGLE::FixGLE(LAMMPS *lmp, int narg, char **arg) : FixGLE::~FixGLE() { delete random; - delete [] A; - delete [] C; - delete [] S; - delete [] T; - delete [] TT; - delete [] ST; + delete[] A; + delete[] C; + delete[] S; + delete[] T; + delete[] TT; + delete[] ST; memory->destroy(sqrt_m); memory->destroy(gle_s); From 0595e4d7cd1f676975b96b71f9fa9f36d904cb09 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 00:24:01 -0500 Subject: [PATCH 6/8] refactor parsing of pair style list parameter file --- src/MISC/pair_list.cpp | 258 +++++++++++++++++------------------------ src/MISC/pair_list.h | 20 ++-- 2 files changed, 113 insertions(+), 165 deletions(-) diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index d1051ecb11..54afcacaed 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -18,19 +18,28 @@ #include "pair_list.h" -#include -#include #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "memory.h" -#include "error.h" +#include "text_file_reader.h" +#include "tokenizer.h" +#include +#include +#include +#include using namespace LAMMPS_NS; -static const char * const stylename[] = { - "none", "harmonic", "morse", "lj126", nullptr +enum { NONE = 0, HARM, MORSE, LJ126 }; + +static std::map stylename = { + { "none", NONE }, + { "harmonic", HARM }, + { "morse", MORSE }, + { "lj126", LJ126 } }; // fast power function for integer exponent > 0 @@ -55,7 +64,6 @@ PairList::PairList(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; respa_enable = 0; cut_global = 0.0; - style = nullptr; params = nullptr; check_flag = 1; } @@ -66,7 +74,6 @@ PairList::~PairList() { memory->destroy(setflag); memory->destroy(cutsq); - memory->destroy(style); memory->destroy(params); } @@ -90,7 +97,7 @@ void PairList::compute(int eflag, int vflag) int pc = 0; for (int n=0; n < npairs; ++n) { - const list_parm_t &par = params[n]; + const list_param &par = params[n]; i = atom->map(par.id1); j = atom->map(par.id2); @@ -122,34 +129,34 @@ void PairList::compute(int eflag, int vflag) if (rsq < par.cutsq) { const double r2inv = 1.0/rsq; - if (style[n] == HARM) { + if (par.style == HARM) { const double r = sqrt(rsq); - const double dr = par.parm.harm.r0 - r; - fpair = 2.0*par.parm.harm.k*dr/r; + const double dr = par.param.harm.r0 - r; + fpair = 2.0*par.param.harm.k*dr/r; if (eflag_either) - epair = par.parm.harm.k*dr*dr - par.offset; + epair = par.param.harm.k*dr*dr - par.offset; - } else if (style[n] == MORSE) { + } else if (par.style == MORSE) { const double r = sqrt(rsq); - const double dr = par.parm.morse.r0 - r; - const double dexp = exp(par.parm.morse.alpha * dr); - fpair = 2.0*par.parm.morse.d0*par.parm.morse.alpha + const double dr = par.param.morse.r0 - r; + const double dexp = exp(par.param.morse.alpha * dr); + fpair = 2.0*par.param.morse.d0*par.param.morse.alpha * (dexp*dexp - dexp) / r; if (eflag_either) - epair = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; + epair = par.param.morse.d0 * (dexp*dexp - 2.0*dexp) - par.offset; - } else if (style[n] == LJ126) { + } else if (par.style == LJ126) { const double r6inv = r2inv*r2inv*r2inv; - const double sig6 = mypow(par.parm.lj126.sigma,6); - fpair = 24.0*par.parm.lj126.epsilon*r6inv + const double sig6 = mypow(par.param.lj126.sigma,6); + fpair = 24.0*par.param.lj126.epsilon*r6inv * (2.0*sig6*sig6*r6inv - sig6) * r2inv; if (eflag_either) - epair = 4.0*par.parm.lj126.epsilon*r6inv + epair = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6) - par.offset; } @@ -210,130 +217,73 @@ void PairList::settings(int narg, char **arg) if (strcmp(arg[2],"check") == 0) check_flag = 1; } - FILE *fp = utils::open_potential(arg[0],lmp,nullptr); - char line[1024]; - if (fp == nullptr) - error->all(FLERR,"Cannot open pair list file"); + std::vector mystyles; + std::vector myparams; - // count lines in file for upper limit of storage needed - int num = 1; - while (fgets(line,1024,fp)) ++num; - rewind(fp); - memory->create(style,num,"pair_list:style"); - memory->create(params,num,"pair_list:params"); - - // now read and parse pair list file for real - npairs = 0; - char *ptr; - tagint id1, id2; - int nharm=0, nmorse=0, nlj126=0; - - while (fgets(line,1024,fp)) { - ptr = strtok(line," \t\n\r\f"); - - // skip empty lines - if (!ptr) continue; - - // skip comment lines starting with # - if (*ptr == '#') continue; - - // get atom ids of pair - id1 = ATOTAGINT(ptr); - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted pair list file"); - id2 = ATOTAGINT(ptr); - - // get potential type - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted pair list file"); - - style[npairs] = NONE; - list_parm_t &par = params[npairs]; - par.id1 = id1; - par.id2 = id2; - - // harmonic potential - if (strcmp(ptr,stylename[HARM]) == 0) { - style[npairs] = HARM; - - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr == nullptr) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); - par.parm.harm.k = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr == nullptr) || (*ptr == '#')) - error->all(FLERR,"Incorrectly formatted harmonic pair parameters"); - par.parm.harm.r0 = utils::numeric(FLERR,ptr,false,lmp); - - ++nharm; - - // morse potential - } else if (strcmp(ptr,stylename[MORSE]) == 0) { - style[npairs] = MORSE; - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.d0 = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.alpha = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted morse pair parameters"); - par.parm.morse.r0 = utils::numeric(FLERR,ptr,false,lmp); - - ++nmorse; - - } else if (strcmp(ptr,stylename[LJ126]) == 0) { - // 12-6 lj potential - style[npairs] = LJ126; - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); - par.parm.lj126.epsilon = utils::numeric(FLERR,ptr,false,lmp); - - ptr = strtok(nullptr," \t\n\r\f"); - if (!ptr) - error->all(FLERR,"Incorrectly formatted 12-6 LJ pair parameters"); - par.parm.lj126.sigma = utils::numeric(FLERR,ptr,false,lmp); - - ++nlj126; - - } else { - error->all(FLERR,"Unknown pair list potential style"); - } - - // optional cutoff parameter. if not specified use global value - ptr = strtok(nullptr," \t\n\r\f"); - if ((ptr != nullptr) && (*ptr != '#')) { - double cut = utils::numeric(FLERR,ptr,false,lmp); - par.cutsq = cut*cut; - } else { - par.cutsq = cut_global*cut_global; - } - - // count complete entry - ++npairs; - } - fclose(fp); - - // informative output + // read and parse potential file only on MPI rank 0. if (comm->me == 0) { - if (screen) - fprintf(screen,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); - if (logfile) - fprintf(logfile,"Read %d (%d/%d/%d) interacting pairs from %s\n", - npairs, nharm, nmorse, nlj126, arg[0]); + int nharm, nmorse, nlj126, nskipped; + FILE *fp = utils::open_potential(arg[0],lmp,nullptr); + TextFileReader reader(fp,"pair list coeffs"); + npairs = nharm = nmorse = nlj126 = nskipped = 0; + + try { + char *line; + while ((line = reader.next_line())) { + ValueTokenizer values(line); + list_param oneparam; + oneparam.id1 = values.next_tagint(); + oneparam.id2 = values.next_tagint(); + oneparam.style = stylename[values.next_string()]; + ++npairs; + + switch (oneparam.style) { + + case HARM: + oneparam.param.harm.k = values.next_double(); + oneparam.param.harm.r0 = values.next_double(); + ++nharm; + break; + + case MORSE: + oneparam.param.morse.d0 = values.next_double(); + oneparam.param.morse.alpha = values.next_double(); + oneparam.param.morse.r0 = values.next_double(); + ++nmorse; + break; + + case LJ126: + oneparam.param.lj126.epsilon = values.next_double(); + oneparam.param.lj126.sigma = values.next_double(); + ++nlj126; + break; + + case NONE: // fallthrough + error->warning(FLERR,"Skipping unrecognized pair list potential entry: {}", + utils::trim(line)); + ++nskipped; + break; + } + if (values.has_next()) + oneparam.cutsq = values.next_double(); + else + oneparam.cutsq = cut_global*cut_global; + + myparams.push_back(oneparam); + } + } catch (std::exception &e) { + error->one(FLERR,"Error reading pair list coeffs file: {}", e.what()); + } + utils::logmesg(lmp, "Read {} ({}/{}/{}) interacting pair lines from {}. " + "{} skipped entries.\n", npairs, nharm, nmorse, nlj126, arg[0], nskipped); + + memory->create(params,npairs,"pair_list:params"); + memcpy(params, myparams.data(),npairs*sizeof(list_param)); + fclose(fp); } + MPI_Bcast(&npairs, 1, MPI_INT, 0, world); + if (comm->me != 0) memory->create(params,npairs,"pair_list:params"); + MPI_Bcast(params, npairs*sizeof(list_param), MPI_CHAR, 0, world); } /* ---------------------------------------------------------------------- @@ -374,21 +324,21 @@ void PairList::init_style() if (offset_flag) { for (int n=0; n < npairs; ++n) { - list_parm_t &par = params[n]; + list_param &par = params[n]; - if (style[n] == HARM) { - const double dr = sqrt(par.cutsq) - par.parm.harm.r0; - par.offset = par.parm.harm.k*dr*dr; + if (par.style == HARM) { + const double dr = sqrt(par.cutsq) - par.param.harm.r0; + par.offset = par.param.harm.k*dr*dr; - } else if (style[n] == MORSE) { - const double dr = par.parm.morse.r0 - sqrt(par.cutsq); - const double dexp = exp(par.parm.morse.alpha * dr); - par.offset = par.parm.morse.d0 * (dexp*dexp - 2.0*dexp); + } else if (par.style == MORSE) { + const double dr = par.param.morse.r0 - sqrt(par.cutsq); + const double dexp = exp(par.param.morse.alpha * dr); + par.offset = par.param.morse.d0 * (dexp*dexp - 2.0*dexp); - } else if (style[n] == LJ126) { + } else if (par.style == LJ126) { const double r6inv = par.cutsq*par.cutsq*par.cutsq; - const double sig6 = mypow(par.parm.lj126.sigma,6); - par.offset = 4.0*par.parm.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); + const double sig6 = mypow(par.param.lj126.sigma,6); + par.offset = 4.0*par.param.lj126.epsilon*r6inv * (sig6*sig6*r6inv - sig6); } } } @@ -411,7 +361,7 @@ double PairList::init_one(int, int) double PairList::memory_usage() { double bytes = (double)npairs * sizeof(int); - bytes += (double)npairs * sizeof(list_parm_t); + bytes += (double)npairs * sizeof(list_param); const int n = atom->ntypes+1; bytes += (double)n*(n*sizeof(int) + sizeof(int *)); bytes += (double)n*(n*sizeof(double) + sizeof(double *)); diff --git a/src/MISC/pair_list.h b/src/MISC/pair_list.h index 60fd3a7a24..b7161c4f2c 100644 --- a/src/MISC/pair_list.h +++ b/src/MISC/pair_list.h @@ -39,8 +39,6 @@ class PairList : public Pair { protected: void allocate(); - enum { NONE = 0, HARM, MORSE, LJ126 }; - // potential specific parameters struct harm_p { double k, r0; @@ -52,23 +50,23 @@ class PairList : public Pair { double epsilon, sigma; }; - union parm_u { - struct harm_p harm; - struct morse_p morse; - struct lj126_p lj126; + union param_u { + harm_p harm; + morse_p morse; + lj126_p lj126; }; - typedef struct { + struct list_param { + int style; // potential style indicator tagint id1, id2; // global atom ids double cutsq; // cutoff**2 for this pair double offset; // energy offset - union parm_u parm; // parameters for style - } list_parm_t; + union param_u param; // parameters for style + }; protected: double cut_global; // global cutoff distance - int *style; // list of styles for pair interactions - list_parm_t *params; // lisf of pair interaction parameters + list_param *params; // lisf of pair interaction parameters int npairs; // # of atom pairs in global list int check_flag; // 1 if checking for missing pairs }; From f5ad91f9feef5b8179cc5e1dd2e8a12c78579102 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 06:28:53 -0500 Subject: [PATCH 7/8] minor tweaks --- src/MANYBODY/pair_tersoff.cpp | 3 +-- src/MISC/pair_list.cpp | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index f21c3bffca..484adc7436 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -431,8 +431,7 @@ void PairTersoff::read_file(char *file) // transparently convert units for supported conversions 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))) { try { ValueTokenizer values(line); diff --git a/src/MISC/pair_list.cpp b/src/MISC/pair_list.cpp index 54afcacaed..5215cdc7e2 100644 --- a/src/MISC/pair_list.cpp +++ b/src/MISC/pair_list.cpp @@ -283,7 +283,7 @@ void PairList::settings(int narg, char **arg) } MPI_Bcast(&npairs, 1, MPI_INT, 0, world); if (comm->me != 0) memory->create(params,npairs,"pair_list:params"); - MPI_Bcast(params, npairs*sizeof(list_param), MPI_CHAR, 0, world); + MPI_Bcast(params, npairs*sizeof(list_param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- From 7676f66674f22296e711a6e358aa93ce424fa09b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 19 Jan 2022 06:29:27 -0500 Subject: [PATCH 8/8] modernize potential file parsing --- src/MANYBODY/pair_edip.cpp | 203 ++++++--------- src/MANYBODY/pair_edip.h | 21 +- src/MANYBODY/pair_edip_multi.cpp | 424 ++++++++++++++----------------- src/MANYBODY/pair_edip_multi.h | 16 +- 4 files changed, 294 insertions(+), 370 deletions(-) diff --git a/src/MANYBODY/pair_edip.cpp b/src/MANYBODY/pair_edip.cpp index af2314b8b9..611b1af455 100644 --- a/src/MANYBODY/pair_edip.cpp +++ b/src/MANYBODY/pair_edip.cpp @@ -31,6 +31,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -368,7 +369,7 @@ void PairEDIP::compute(int eflag, int vflag) directorCos_ik_z = invR_ik * dr_ik[2]; cosTeta = directorCos_ij_x * directorCos_ik_x + directorCos_ij_y * directorCos_ik_y + - directorCos_ij_z * directorCos_ik_z; + directorCos_ij_z * directorCos_ik_z; cosTetaDiff = cosTeta + tauFunction; cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff; @@ -376,33 +377,33 @@ void PairEDIP::compute(int eflag, int vflag) expMinusQFunctionCosTetaDiffCosTetaDiff = exp(-qFunctionCosTetaDiffCosTetaDiff); potentia3B_factor = lambda * - ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + - eta * qFunctionCosTetaDiffCosTetaDiff); + ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) + + eta * qFunctionCosTetaDiffCosTetaDiff); exp3B_ik = preExp3B_ij[neighbor_k]; exp3BDerived_ik = preExp3BDerived_ij[neighbor_k]; forceMod3B_factor1_ij = -exp3BDerived_ij * exp3B_ik * potentia3B_factor; forceMod3B_factor2 = 2.0 * lambda * exp3B_ij * exp3B_ik * qFunction * cosTetaDiff * - (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); + (eta + expMinusQFunctionCosTetaDiffCosTetaDiff); forceMod3B_factor2_ij = forceMod3B_factor2 * invR_ij; f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_x - directorCos_ik_x); f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_y - directorCos_ik_y); f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z + - forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); + forceMod3B_factor2_ij * (cosTeta * directorCos_ij_z - directorCos_ik_z); forceMod3B_factor1_ik = -exp3BDerived_ik * exp3B_ij * potentia3B_factor; forceMod3B_factor2_ik = forceMod3B_factor2 * invR_ik; f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_x - directorCos_ij_x); f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_y - directorCos_ij_y); f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z + - forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); + forceMod3B_factor2_ik * (cosTeta * directorCos_ik_z - directorCos_ij_z); forceModCoord += (forceMod3B_factor2 * (tauFunctionDerived - 0.5 * mu * cosTetaDiff)); @@ -765,136 +766,92 @@ double PairEDIP::init_one(int i, int j) void PairEDIP::read_file(char *file) { - int params_per_line = 20; - char **words = new char *[params_per_line + 1]; - memory->sfree(params); params = nullptr; nparams = maxparam = 0; // open file on proc 0 - FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file, lmp, nullptr); - if (fp == nullptr) - error->one(FLERR, "Cannot open EDIP potential file {}: {}", file, utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "edip"); + char *line; - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + while ((line = reader.next_line(NPARAMS_PER_LINE))) { + try { + ValueTokenizer values(line); + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - int n, nwords, ielement, jelement, kelement; - char line[MAXLINE], *ptr; - int eof = 0; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - while (true) { - if (comm->me == 0) { - 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); + 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; - // strip comment, skip line if blank + // load up parameter settings and error check their values - if ((ptr = strchr(line, '#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); - // concatenate additional lines until have params_per_line words + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n], MAXLINE - n, fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else - n = strlen(line) + 1; + memset(params + nparams, 0, DELTA*sizeof(Param)); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].A = values.next_double(); + params[nparams].B = values.next_double(); + params[nparams].cutoffA = values.next_double(); + params[nparams].cutoffC = values.next_double(); + params[nparams].alpha = values.next_double(); + params[nparams].beta = values.next_double(); + params[nparams].eta = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lambda = values.next_double(); + params[nparams].mu = values.next_double(); + params[nparams].rho = values.next_double(); + params[nparams].sigma = values.next_double(); + params[nparams].Q0 = values.next_double(); + params[nparams].u1 = values.next_double(); + params[nparams].u2 = values.next_double(); + params[nparams].u3 = values.next_double(); + params[nparams].u4 = values.next_double(); + } catch (std::exception &e) { + error->one(FLERR, "Error reading EDIP potential file: {}", e.what()); } - 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); - if ((ptr = strchr(line, '#'))) *ptr = '\0'; - nwords = utils::count_words(line); + + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); + + nparams++; } - - if (nwords != params_per_line) error->all(FLERR, "Incorrect format in EDIP potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line, " \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr, " \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 entry in file - - 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"); - - // 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; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].A = atof(words[3]); - params[nparams].B = atof(words[4]); - params[nparams].cutoffA = atof(words[5]); - params[nparams].cutoffC = atof(words[6]); - params[nparams].alpha = atof(words[7]); - params[nparams].beta = atof(words[8]); - params[nparams].eta = atof(words[9]); - params[nparams].gamm = atof(words[10]); - params[nparams].lambda = atof(words[11]); - params[nparams].mu = atof(words[12]); - params[nparams].rho = atof(words[13]); - params[nparams].sigma = atof(words[14]); - params[nparams].Q0 = atof(words[15]); - params[nparams].u1 = atof(words[16]); - params[nparams].u2 = atof(words[17]); - params[nparams].u3 = atof(words[18]); - params[nparams].u4 = atof(words[19]); - - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || - params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || - params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamm < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || - params[nparams].sigma < 0.0) - error->all(FLERR, "Illegal EDIP parameter"); - - nparams++; } + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete[] words; + if (comm->me != 0) + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -946,7 +903,7 @@ void PairEDIP::setup_params() cutoffC = params[0].cutoffC; sigma = params[0].sigma; lambda = params[0].lambda; - gamm = params[0].gamm; + gamm = params[0].gamma; eta = params[0].eta; Q0 = params[0].Q0; mu = params[0].mu; diff --git a/src/MANYBODY/pair_edip.h b/src/MANYBODY/pair_edip.h index 5812768d55..a0d1e45613 100644 --- a/src/MANYBODY/pair_edip.h +++ b/src/MANYBODY/pair_edip.h @@ -34,14 +34,23 @@ class PairEDIP : public Pair { double init_one(int, int); void init_style(); + static constexpr int NPARAMS_PER_LINE = 20; + protected: struct Param { - double A, B; - double cutoffA, cutoffC, cutsq; - double alpha, beta; - double eta, gamm, lambda, mu, rho, sigma, Q0; - double u1, u2, u3, u4; - int ielement, jelement, kelement; + double A, B; // coefficients for pair interaction I-J + double cutoffA; // cut-off distance for pair interaction I-J + double cutoffC; // lower cut-off distance for calculating Z_I + double alpha; // coefficient for calculating Z_I + double beta; // attractive term for pair I-J + double sigma; // cut-off coefficient for pair I-J + double rho; // pair I-J + double gamma; // coefficient for three-body interaction I-J-K + double eta, lambda; // coefficients for function h(l,Z) + double mu, Q0; // coefficients for function Q(Z) + double u1, u2, u3, u4; // coefficients for function tau(Z) + double cutsq; + int ielement, jelement, kelement, dummy; // dummy added for better alignment }; double *preInvR_ij; diff --git a/src/MANYBODY/pair_edip_multi.cpp b/src/MANYBODY/pair_edip_multi.cpp index 8017fa4f8e..ffbc135a34 100644 --- a/src/MANYBODY/pair_edip_multi.cpp +++ b/src/MANYBODY/pair_edip_multi.cpp @@ -30,6 +30,7 @@ #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -118,8 +119,8 @@ void PairEDIPMulti::compute(int eflag, int vflag) double costheta; double dpairZ,dtripleZ; - // eflag != 0 means compute energy contributions in this step - // vflag != 0 means compute virial contributions in this step + // eflag != 0 means compute energy contributions in this step + // vflag != 0 means compute virial contributions in this step evdwl = 0.0; ev_init(eflag,vflag); @@ -155,39 +156,40 @@ void PairEDIPMulti::compute(int eflag, int vflag) // pre-loop to compute environment coordination f(Z) for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; + j = jlist[jj]; + j &= NEIGHMASK; - double delx, dely, delz, r_ij; + double delx, dely, delz, r_ij; - delx = x[j][0] - xtmp; - dely = x[j][1] - ytmp; - delz = x[j][2] - ztmp; - r_ij = delx * delx + dely * dely + delz * delz; + delx = x[j][0] - xtmp; + dely = x[j][1] - ytmp; + delz = x[j][2] - ztmp; + r_ij = delx * delx + dely * dely + delz * delz; - jtype = map[type[j]]; - ijparam = elem3param[itype][jtype][jtype]; - if (r_ij > params[ijparam].cutsq) continue; + jtype = map[type[j]]; + const Param ¶m = params[elem3param[itype][jtype][jtype]]; - r_ij = sqrt(r_ij); + if (r_ij > param.cutsq) continue; - // zeta and its derivative dZ/dr + r_ij = sqrt(r_ij); - if (r_ij < params[ijparam].cutoffC) zeta_i += 1.0; - else { - double f, fdr; - edip_fc(r_ij, ¶ms[ijparam], f, fdr); - zeta_i += f; - dzetair = -fdr / r_ij; + // zeta and its derivative dZ/dr - preForceCoord_counter=numForceCoordPairs*5; - preForceCoord[preForceCoord_counter+0]=dzetair; - preForceCoord[preForceCoord_counter+1]=delx; - preForceCoord[preForceCoord_counter+2]=dely; - preForceCoord[preForceCoord_counter+3]=delz; - preForceCoord[preForceCoord_counter+4]=j; - numForceCoordPairs++; - } + if (r_ij < param.cutoffC) zeta_i += 1.0; + else { + double f, fdr; + edip_fc(r_ij, param, f, fdr); + zeta_i += f; + dzetair = -fdr / r_ij; + + preForceCoord_counter=numForceCoordPairs*5; + preForceCoord[preForceCoord_counter+0]=dzetair; + preForceCoord[preForceCoord_counter+1]=delx; + preForceCoord[preForceCoord_counter+2]=dely; + preForceCoord[preForceCoord_counter+3]=delz; + preForceCoord[preForceCoord_counter+4]=j; + numForceCoordPairs++; + } } // two-body interactions @@ -217,7 +219,7 @@ void PairEDIPMulti::compute(int eflag, int vflag) // already considered in constructing the potential double fdr, fdZ; - edip_pair(r_ij, zeta_i, ¶ms[ijparam], evdwl, fdr, fdZ); + edip_pair(r_ij, zeta_i, params[ijparam], evdwl, fdr, fdZ); fpair = -fdr / r_ij; dpairZ += fdZ; @@ -234,102 +236,102 @@ void PairEDIPMulti::compute(int eflag, int vflag) // three-body Forces for (kk = jj + 1; kk < jnum; kk++) { - double dr_ik[3], r_ik, f_ik[3]; + double dr_ik[3], r_ik, f_ik[3]; - k = jlist[kk]; - k &= NEIGHMASK; - ktype = map[type[k]]; - ikparam = elem3param[itype][ktype][ktype]; - ijkparam = elem3param[itype][jtype][ktype]; + k = jlist[kk]; + k &= NEIGHMASK; + ktype = map[type[k]]; + ikparam = elem3param[itype][ktype][ktype]; + ijkparam = elem3param[itype][jtype][ktype]; - dr_ik[0] = x[k][0] - xtmp; - dr_ik[1] = x[k][1] - ytmp; - dr_ik[2] = x[k][2] - ztmp; - r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; + dr_ik[0] = x[k][0] - xtmp; + dr_ik[1] = x[k][1] - ytmp; + dr_ik[2] = x[k][2] - ztmp; + r_ik = dr_ik[0]*dr_ik[0] + dr_ik[1]*dr_ik[1] + dr_ik[2]*dr_ik[2]; - if (r_ik > params[ikparam].cutsq) continue; + if (r_ik > params[ikparam].cutsq) continue; - r_ik = sqrt(r_ik); + r_ik = sqrt(r_ik); - costheta=dot3(dr_ij, dr_ik) / r_ij / r_ik; + costheta=dot3(dr_ij, dr_ik) / r_ij / r_ik; - double v1, v2, v3, v4, v5, v6, v7; + double v1, v2, v3, v4, v5, v6, v7; - edip_fcut3(r_ij, ¶ms[ijparam], v1, v2); - edip_fcut3(r_ik, ¶ms[ikparam], v3, v4); - edip_h(costheta, zeta_i, ¶ms[ijkparam], v5, v6, v7); + edip_fcut3(r_ij, params[ijparam], v1, v2); + edip_fcut3(r_ik, params[ikparam], v3, v4); + edip_h(costheta, zeta_i, params[ijkparam], v5, v6, v7); - // potential energy and forces - evdwl = v1 * v3 * v5; - dtripleZ += v1 * v3 * v7; + // potential energy and forces + evdwl = v1 * v3 * v5; + dtripleZ += v1 * v3 * v7; - double dri[3], drj[3], drk[3]; - double dhl, dfr; + double dri[3], drj[3], drk[3]; + double dhl, dfr; - dhl = v1 * v3 * v6; + dhl = v1 * v3 * v6; - costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk); + costheta_d(dr_ij, r_ij, dr_ik, r_ik, dri, drj, drk); - f_ij[0] = -dhl * drj[0]; - f_ij[1] = -dhl * drj[1]; - f_ij[2] = -dhl * drj[2]; - f_ik[0] = -dhl * drk[0]; - f_ik[1] = -dhl * drk[1]; - f_ik[2] = -dhl * drk[2]; + f_ij[0] = -dhl * drj[0]; + f_ij[1] = -dhl * drj[1]; + f_ij[2] = -dhl * drj[2]; + f_ik[0] = -dhl * drk[0]; + f_ik[1] = -dhl * drk[1]; + f_ik[2] = -dhl * drk[2]; - dfr = v2 * v3 * v5; - fpair = -dfr / r_ij; + dfr = v2 * v3 * v5; + fpair = -dfr / r_ij; - f_ij[0] += fpair * dr_ij[0]; - f_ij[1] += fpair * dr_ij[1]; - f_ij[2] += fpair * dr_ij[2]; + f_ij[0] += fpair * dr_ij[0]; + f_ij[1] += fpair * dr_ij[1]; + f_ij[2] += fpair * dr_ij[2]; - dfr = v1 * v4 * v5; - fpair = -dfr / r_ik; + dfr = v1 * v4 * v5; + fpair = -dfr / r_ik; - f_ik[0] += fpair * dr_ik[0]; - f_ik[1] += fpair * dr_ik[1]; - f_ik[2] += fpair * dr_ik[2]; + f_ik[0] += fpair * dr_ik[0]; + f_ik[1] += fpair * dr_ik[1]; + f_ik[2] += fpair * dr_ik[2]; - f[j][0] += f_ij[0]; - f[j][1] += f_ij[1]; - f[j][2] += f_ij[2]; + f[j][0] += f_ij[0]; + f[j][1] += f_ij[1]; + f[j][2] += f_ij[2]; - f[k][0] += f_ik[0]; - f[k][1] += f_ik[1]; - f[k][2] += f_ik[2]; + f[k][0] += f_ik[0]; + f[k][1] += f_ik[1]; + f[k][2] += f_ik[2]; - f[i][0] -= f_ij[0] + f_ik[0]; - f[i][1] -= f_ij[1] + f_ik[1]; - f[i][2] -= f_ij[2] + f_ik[2]; + f[i][0] -= f_ij[0] + f_ik[0]; + f[i][1] -= f_ij[1] + f_ik[1]; + f[i][2] -= f_ij[2] + f_ik[2]; - if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); + if (evflag) ev_tally3(i,j,k,evdwl,0.0,f_ij,f_ik,dr_ij,dr_ik); } } // forces due to environment coordination f(Z) for (int idx = 0; idx < numForceCoordPairs; idx++) { - double delx, dely, delz; + double delx, dely, delz; - preForceCoord_counter = idx * 5; - dzetair = preForceCoord[preForceCoord_counter+0]; - delx = preForceCoord[preForceCoord_counter+1]; - dely = preForceCoord[preForceCoord_counter+2]; - delz = preForceCoord[preForceCoord_counter+3]; - j = static_cast (preForceCoord[preForceCoord_counter+4]); + preForceCoord_counter = idx * 5; + dzetair = preForceCoord[preForceCoord_counter+0]; + delx = preForceCoord[preForceCoord_counter+1]; + dely = preForceCoord[preForceCoord_counter+2]; + delz = preForceCoord[preForceCoord_counter+3]; + j = static_cast (preForceCoord[preForceCoord_counter+4]); - dzetair *= (dpairZ + dtripleZ); + dzetair *= (dpairZ + dtripleZ); - f[j][0] += dzetair * delx; - f[j][1] += dzetair * dely; - f[j][2] += dzetair * delz; + f[j][0] += dzetair * delx; + f[j][1] += dzetair * dely; + f[j][2] += dzetair * delz; - f[i][0] -= dzetair * delx; - f[i][1] -= dzetair * dely; - f[i][2] -= dzetair * delz; + f[i][0] -= dzetair * delx; + f[i][1] -= dzetair * dely; + f[i][2] -= dzetair * delz; - evdwl = 0.0; - if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz); + evdwl = 0.0; + if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, dzetair, -delx, -dely, -delz); } } @@ -342,13 +344,13 @@ double sqr(double x) } //pair Vij, partial derivatives dVij(r,Z)/dr and dVij(r,Z)/dZ -void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng, +void PairEDIPMulti::edip_pair(double r, double z, const Param ¶m, double &eng, double &fdr, double &fZ) { - double A = param->A; - double B = param->B; - double rho = param->rho; - double beta = param->beta; + double A = param.A; + double B = param.B; + double rho = param.rho; + double beta = param.beta; double v1,v2,v3,v4; v1 = pow(B / r, rho); @@ -361,11 +363,11 @@ void PairEDIPMulti::edip_pair(double r, double z, Param *param, double &eng, } //function fc(r) in calculating coordination Z and derivative fc'(r) -void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fc(double r, const Param ¶m, double &f, double &fdr) { - double a = param->cutoffA; - double c = param->cutoffC; - double alpha = param->alpha; + double a = param.cutoffA; + double c = param.cutoffC; + double alpha = param.alpha; double x; double v1, v2; @@ -392,10 +394,10 @@ void PairEDIPMulti::edip_fc(double r, Param *param, double &f, double &fdr) } //cut-off function for Vij and its derivative fcut2'(r) -void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fcut2(double r, const Param ¶m, double &f, double &fdr) { - double sigma = param->sigma; - double a = param->cutoffA; + double sigma = param.sigma; + double a = param.cutoffA; double v1; if (r > a - 1E-6) @@ -411,12 +413,12 @@ void PairEDIPMulti::edip_fcut2(double r, Param *param, double &f, double &fdr) } //function tau(Z) and its derivative tau'(Z) -void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ) +void PairEDIPMulti::edip_tau(double z, const Param ¶m, double &f, double &fdZ) { - double u1 = param->u1; - double u2 = param->u2; - double u3 = param->u3; - double u4 = param->u4; + double u1 = param.u1; + double u2 = param.u2; + double u3 = param.u3; + double u4 = param.u4; double v1, v2; v1 = exp(-u4 * z); @@ -427,13 +429,13 @@ void PairEDIPMulti::edip_tau(double z, Param *param, double &f, double &fdZ) } //function h(l,Z) and its partial derivatives dh(l,Z)/dl and dh(l,Z)/dZ -void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f, +void PairEDIPMulti::edip_h(double l, double z, const Param ¶m, double &f, double &fdl, double &fdZ) { - double lambda = param->lambda; - double eta = param->eta; - double Q0 = param->Q0; - double mu = param->mu; + double lambda = param.lambda; + double eta = param.eta; + double Q0 = param.Q0; + double mu = param.mu; double Q, QdZ, Tau, TaudZ; double u2, du2l, du2Z; double v1, v2, v3; @@ -464,10 +466,10 @@ void PairEDIPMulti::edip_h(double l, double z, Param *param, double &f, } //cut-off function for Vijk and its derivative fcut3'(r) -void PairEDIPMulti::edip_fcut3(double r, Param *param, double &f, double &fdr) +void PairEDIPMulti::edip_fcut3(double r, const Param ¶m, double &f, double &fdr) { - double gamma = param->gamma; - double a = param->cutoffA; + double gamma = param.gamma; + double a = param.cutoffA; double v1; if (r > a - 1E-6) @@ -578,137 +580,92 @@ double PairEDIPMulti::init_one(int i, int j) void PairEDIPMulti::read_file(char *file) { - int params_per_line = 20; - char **words = new char*[params_per_line+1]; - memory->sfree(params); params = nullptr; nparams = maxparam = 0; // open file on proc 0 - FILE *fp; if (comm->me == 0) { - fp = utils::open_potential(file,lmp,nullptr); - if (fp == nullptr) - error->one(FLERR,"Cannot open EDIP potential file {}: {}",file,utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "edip"); + char *line; - // read each set of params from potential file - // one set of params can span multiple lines - // store params if all 3 element tags are in element list + while ((line = reader.next_line(NPARAMS_PER_LINE))) { + try { + ValueTokenizer values(line); + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; + // ielement,jelement,kelement = 1st args + // if all 3 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement, kelement; - while (true) { - if (comm->me == 0) { - 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); + 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; - // strip comment, skip line if blank + // load up parameter settings and error check their values - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); - if (nwords == 0) continue; + if (nparams == maxparam) { + maxparam += DELTA; + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), + "pair:params"); - // concatenate additional lines until have params_per_line words + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind - while (nwords < params_per_line) { - n = strlen(line); - if (comm->me == 0) { - ptr = fgets(&line[n],MAXLINE-n,fp); - if (ptr == nullptr) { - eof = 1; - fclose(fp); - } else n = strlen(line) + 1; + memset(params + nparams, 0, DELTA*sizeof(Param)); + } + + params[nparams].ielement = ielement; + params[nparams].jelement = jelement; + params[nparams].kelement = kelement; + params[nparams].A = values.next_double(); + params[nparams].B = values.next_double(); + params[nparams].cutoffA = values.next_double(); + params[nparams].cutoffC = values.next_double(); + params[nparams].alpha = values.next_double(); + params[nparams].beta = values.next_double(); + params[nparams].eta = values.next_double(); + params[nparams].gamma = values.next_double(); + params[nparams].lambda = values.next_double(); + params[nparams].mu = values.next_double(); + params[nparams].rho = values.next_double(); + params[nparams].sigma = values.next_double(); + params[nparams].Q0 = values.next_double(); + params[nparams].u1 = values.next_double(); + params[nparams].u2 = values.next_double(); + params[nparams].u3 = values.next_double(); + params[nparams].u4 = values.next_double(); + } catch (std::exception &e) { + error->one(FLERR, "Error reading EDIP potential file: {}", e.what()); } - 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); - if ((ptr = strchr(line,'#'))) *ptr = '\0'; - nwords = utils::count_words(line); + + if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || params[nparams].cutoffA < 0.0 || + params[nparams].cutoffC < 0.0 || params[nparams].alpha < 0.0 || + params[nparams].beta < 0.0 || params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || + params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || params[nparams].rho < 0.0 || + params[nparams].sigma < 0.0) + error->all(FLERR, "Illegal EDIP parameter"); + + nparams++; } - - if (nwords != params_per_line) - error->all(FLERR,"Incorrect format in EDIP potential file"); - - // words = ptrs to all words in line - - nwords = 0; - words[nwords++] = strtok(line," \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr," \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 entry in file - - 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"); - - // 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; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].A = atof(words[3]); - params[nparams].B = atof(words[4]); - params[nparams].cutoffA = atof(words[5]); - params[nparams].cutoffC = atof(words[6]); - params[nparams].alpha = atof(words[7]); - params[nparams].beta = atof(words[8]); - params[nparams].eta = atof(words[9]); - params[nparams].gamma = atof(words[10]); - params[nparams].lambda = atof(words[11]); - params[nparams].mu = atof(words[12]); - params[nparams].rho = atof(words[13]); - params[nparams].sigma = atof(words[14]); - params[nparams].Q0 = atof(words[15]); - params[nparams].u1 = atof(words[16]); - params[nparams].u2 = atof(words[17]); - params[nparams].u3 = atof(words[18]); - params[nparams].u4 = atof(words[19]); - - if (params[nparams].A < 0.0 || params[nparams].B < 0.0 || - params[nparams].cutoffA < 0.0 || params[nparams].cutoffC < 0.0 || - params[nparams].alpha < 0.0 || params[nparams].beta < 0.0 || - params[nparams].eta < 0.0 || params[nparams].gamma < 0.0 || - params[nparams].lambda < 0.0 || params[nparams].mu < 0.0 || - params[nparams].rho < 0.0 || params[nparams].sigma < 0.0) - error->all(FLERR,"Illegal EDIP parameter"); - - nparams++; } + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete [] words; + if (comm->me != 0) + params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); + + MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -753,5 +710,4 @@ void PairEDIPMulti::setup() rtmp = sqrt(params[m].cutsq); if (rtmp > cutmax) cutmax = rtmp; } - } diff --git a/src/MANYBODY/pair_edip_multi.h b/src/MANYBODY/pair_edip_multi.h index 3ee7347a56..1070956bc2 100644 --- a/src/MANYBODY/pair_edip_multi.h +++ b/src/MANYBODY/pair_edip_multi.h @@ -34,6 +34,8 @@ class PairEDIPMulti : public Pair { double init_one(int, int); void init_style(); + static constexpr int NPARAMS_PER_LINE = 20; + protected: struct Param { double A, B; // coefficients for pair interaction I-J @@ -48,7 +50,7 @@ class PairEDIPMulti : public Pair { double mu, Q0; // coefficients for function Q(Z) double u1, u2, u3, u4; // coefficients for function tau(Z) double cutsq; - int ielement, jelement, kelement; + int ielement, jelement, kelement, dummy; // dummy added for better alignment }; double *preForceCoord; @@ -63,12 +65,12 @@ class PairEDIPMulti : public Pair { void read_file(char *); void setup(); - void edip_pair(double, double, Param *, double &, double &, double &); - void edip_fc(double, Param *, double &, double &); - void edip_fcut2(double, Param *, double &, double &); - void edip_tau(double, Param *, double &, double &); - void edip_h(double, double, Param *, double &, double &, double &); - void edip_fcut3(double, Param *, double &, double &); + void edip_pair(double, double, const Param &, double &, double &, double &); + void edip_fc(double, const Param &, double &, double &); + void edip_fcut2(double, const Param &, double &, double &); + void edip_tau(double, const Param &, double &, double &); + void edip_h(double, double, const Param &, double &, double &, double &); + void edip_fcut3(double, const Param &, double &, double &); }; } // namespace LAMMPS_NS