Add PotentialFileReader and use it in PairSW

This commit is contained in:
Richard Berger 2020-05-26 16:45:43 -04:00
parent 74249380ec
commit c302c0bca2
No known key found for this signature in database
GPG Key ID: A9E83994E0BA0CAB
4 changed files with 187 additions and 104 deletions

View File

@ -29,6 +29,8 @@
#include "neigh_list.h"
#include "neigh_request.h"
#include "utils.h"
#include "tokenizer.h"
#include "potential_file_reader.h"
using namespace LAMMPS_NS;
@ -347,128 +349,84 @@ double PairSW::init_one(int i, int j)
void PairSW::read_file(char *file)
{
int params_per_line = 14;
char **words = new char*[params_per_line+1];
memory->sfree(params);
params = NULL;
nparams = maxparam = 0;
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
snprintf(str,128,"Cannot open Stillinger-Weber potential file %s",file);
error->one(FLERR,str);
}
}
PotentialFileReader reader(lmp, file, "Stillinger-Weber");
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(Param::NPARAMS_PER_LINE)) {
try {
ValueTokenizer values(line, " \t\n\r\f");
int n,nwords,ielement,jelement,kelement;
char line[MAXLINE],*ptr;
int eof = 0;
std::string iword = values.next_string();
std::string jword = values.next_string();
std::string kword = values.next_string();
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
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);
// 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;
// strip comment, skip line if blank
for (ielement = 0; ielement < nelements; ielement++)
if (iword == elements[ielement]) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (jword == elements[jelement]) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (kword == elements[kelement]) break;
if (kelement == nelements) continue;
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = utils::count_words(line);
if (nwords == 0) continue;
// load up parameter settings and error check their values
// concatenate additional lines until have params_per_line words
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].epsilon = values.next_double();
params[nparams].sigma = values.next_double();
params[nparams].littlea = values.next_double();
params[nparams].lambda = values.next_double();
params[nparams].gamma = values.next_double();
params[nparams].costheta = values.next_double();
params[nparams].biga = values.next_double();
params[nparams].bigb = values.next_double();
params[nparams].powerp = values.next_double();
params[nparams].powerq = values.next_double();
params[nparams].tol = values.next_double();
} 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].epsilon < 0.0 || params[nparams].sigma < 0.0 ||
params[nparams].littlea < 0.0 || params[nparams].lambda < 0.0 ||
params[nparams].gamma < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].bigb < 0.0 || params[nparams].powerp < 0.0 ||
params[nparams].powerq < 0.0 || params[nparams].tol < 0.0)
error->one(FLERR,"Illegal Stillinger-Weber parameter");
nparams++;
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in Stillinger-Weber potential file");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \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");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].epsilon = atof(words[3]);
params[nparams].sigma = atof(words[4]);
params[nparams].littlea = atof(words[5]);
params[nparams].lambda = atof(words[6]);
params[nparams].gamma = atof(words[7]);
params[nparams].costheta = atof(words[8]);
params[nparams].biga = atof(words[9]);
params[nparams].bigb = atof(words[10]);
params[nparams].powerp = atof(words[11]);
params[nparams].powerq = atof(words[12]);
params[nparams].tol = atof(words[13]);
if (params[nparams].epsilon < 0.0 || params[nparams].sigma < 0.0 ||
params[nparams].littlea < 0.0 || params[nparams].lambda < 0.0 ||
params[nparams].gamma < 0.0 || params[nparams].biga < 0.0 ||
params[nparams].bigb < 0.0 || params[nparams].powerp < 0.0 ||
params[nparams].powerq < 0.0 || params[nparams].tol < 0.0)
error->all(FLERR,"Illegal Stillinger-Weber parameter");
nparams++;
}
delete [] words;
MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
if(comm->me != 0) {
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
}
MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
}
/* ---------------------------------------------------------------------- */

View File

@ -44,6 +44,7 @@ class PairSW : public Pair {
double sigma_gamma,lambda_epsilon,lambda_epsilon2;
double c1,c2,c3,c4,c5,c6;
int ielement,jelement,kelement;
static const int NPARAMS_PER_LINE = 14;
};
protected:

View File

@ -0,0 +1,82 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Richard Berger (Temple U)
------------------------------------------------------------------------- */
#include "lammps.h"
#include "force.h"
#include "error.h"
#include "comm.h"
#include "potential_file_reader.h"
#include "utils.h"
#include <cstring>
using namespace LAMMPS_NS;
PotentialFileReader::PotentialFileReader(LAMMPS *lmp,
const std::string &filename,
const std::string &potential_name) : Pointers(lmp) {
if (comm->me != 0) {
error->one(FLERR, "PotentialFileReader should only be called by proc 0!");
}
fp = force->open_potential(filename.c_str());
if (fp == NULL) {
char str[128];
snprintf(str, 128, "cannot open %s potential file %s", potential_name.c_str(), filename.c_str());
error->one(FLERR, str);
}
}
PotentialFileReader::~PotentialFileReader() {
fclose(fp);
}
char *PotentialFileReader::next_line(int nparams) {
// concatenate lines until have nparams words
int n = 0;
int nwords = 0;
do {
char *ptr = fgets(&line[n], MAXLINE - n, fp);
if (ptr == nullptr) {
// EOF
if (nwords > 0 && nwords < nparams)
{
char str[128];
snprintf(str, 128, "Incorrect format in %s potential file", potential_name.c_str());
error->all(FLERR, str);
}
return nullptr;
}
// strip comment
if ((ptr = strchr(line, '#'))) *ptr = '\0';
nwords = utils::count_words(line);
// skip line if blank
if (nwords == 0)
continue;
n = strlen(line) + 1;
} while (nwords < nparams);
return line;
}

View File

@ -0,0 +1,42 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Richard Berger (Temple U)
------------------------------------------------------------------------- */
#ifndef LMP_POTENTIAL_FILE_READER_H
#define LMP_POTENTIAL_FILE_READER_H
#include <string>
#include "pointers.h"
namespace LAMMPS_NS
{
class PotentialFileReader : protected Pointers {
std::string potential_name;
static const int MAXLINE = 1024;
char line[MAXLINE];
FILE *fp;
public:
PotentialFileReader(class LAMMPS *lmp, const std::string &filename, const std::string &potential_name);
~PotentialFileReader();
char *next_line(int nparams);
};
} // namespace LAMMPS_NS
#endif