diff --git a/examples/PACKAGES/extep/BN.extep b/examples/PACKAGES/extep/BN.extep new file mode 120000 index 0000000000..2f52a7ec7b --- /dev/null +++ b/examples/PACKAGES/extep/BN.extep @@ -0,0 +1 @@ +../../../potentials/BN.extep \ No newline at end of file diff --git a/potentials/BN.extep b/potentials/BN.extep index 8732ada84b..967a0895bf 100644 --- a/potentials/BN.extep +++ b/potentials/BN.extep @@ -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 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 -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 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 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 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 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.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.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.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 # diff --git a/src/MANYBODY/pair_extep.cpp b/src/MANYBODY/pair_extep.cpp index 7c75cfcfc0..4e7de1dadd 100644 --- a/src/MANYBODY/pair_extep.cpp +++ b/src/MANYBODY/pair_extep.cpp @@ -28,6 +28,7 @@ #include "my_page.h" #include "neigh_list.h" #include "neighbor.h" +#include "potential_file_reader.h" #include #include @@ -511,234 +512,177 @@ double PairExTeP::init_one(int i, int j) void PairExTeP::read_file(char *file) { - int params_per_line = 17; - 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 ExTeP potential file {}: {}",file,utils::getsyserror()); - } + PotentialFileReader reader(lmp, file, "ExTeP", unit_convert_flag); + char *line; - // read each line out of file, skipping blank lines or leading '#' - // store line of params if all 3 element tags are in element list + // transparently convert units for supported conversions - int n,nwords,ielement,jelement,kelement; - char line[MAXLINE],*ptr; - int eof = 0; + int unit_convert = reader.get_unit_convert(); + double conversion_factor = utils::get_conversion_factor(utils::ENERGY,unit_convert); - 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); + while ((line = reader.next_line(17))) { + try { + ValueTokenizer values(line); - // 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'; - nwords = utils::count_words(line); - if (nwords == 0) 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 + 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) { - 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; + // 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].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; - 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].c < 0.0 || + params[nparams].d < 0.0 || + params[nparams].powern < 0.0 || + params[nparams].beta < 0.0 || + 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) - error->all(FLERR,"Insufficient spline parameters in potential file"); + /* F_IJ (3) */ + // 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; - words[nwords++] = strtok(line," \t\n\r\f"); - while ((words[nwords++] = strtok(nullptr," \t\n\r\f"))) continue; + while ((line = reader.next_line(8))) { + try { + ValueTokenizer values(line); - // ielement,jelement,kelement = 1st args - // if all 3 args are in element list, then parse this line - // else skip to next line + std::string iname = values.next_string(); + std::string jname = values.next_string(); + std::string kname = values.next_string(); - 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; + // ielement,jelement = 1st args + // if all 2 args are in element list, then parse this line + // else skip to next entry in file + int ielement, jelement; - // 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) { - maxparam += DELTA; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); + // 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. - // make certain all addional allocated storage is initialized - // to avoid false positives when checking with valgrind + if (!utils::is_integer(kname)) + 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; - params[nparams].jelement = jelement; - params[nparams].kelement = kelement; - params[nparams].powerm = atof(words[3]); - params[nparams].gamma = atof(words[4]); - params[nparams].lam3 = atof(words[5]); - 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]); + // 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; - // currently only allow m exponent of 1 or 3 - - params[nparams].powermint = int(params[nparams].powerm); - - if (params[nparams].c < 0.0 || params[nparams].d < 0.0 || - params[nparams].powern < 0.0 || params[nparams].beta < 0.0 || - 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->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;ielme == 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; + 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; + } + } + } catch (TokenizerException &e) { + error->one(FLERR, e.what()); } } } + + MPI_Bcast(&nparams, 1, MPI_INT, 0, world); + MPI_Bcast(&maxparam, 1, MPI_INT, 0, world); - delete[] words; - /* END F_IJ (3) */ - + if (comm->me != 0) { + 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); } /* ---------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index c5578dfc18..37ea0bfebf 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -428,6 +428,7 @@ void PairTersoff::read_file(char *file) int unit_convert = reader.get_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); @@ -455,8 +456,7 @@ void PairTersoff::read_file(char *file) if (nparams == maxparam) { maxparam += DELTA; - params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), - "pair:params"); + 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 @@ -505,7 +505,7 @@ void PairTersoff::read_file(char *file) params[nparams].biga < 0.0 || params[nparams].powerm - params[nparams].powermint != 0.0 || (params[nparams].powermint != 3 && - params[nparams].powermint != 1) || + params[nparams].powermint != 1) || params[nparams].gamma < 0.0) error->one(FLERR,"Illegal Tersoff parameter");