Merge branch 'ml-pod-fixes' into collected-small-changes

This commit is contained in:
Axel Kohlmeyer
2022-12-05 19:18:03 -05:00
3 changed files with 60 additions and 55 deletions

View File

@ -46,6 +46,10 @@ using MathSpecial::powint;
static constexpr double SMALL = 1.0e-10;
FitPOD::FitPOD(LAMMPS *_lmp) : Command(_lmp), podptr(nullptr)
{
}
void FitPOD::command(int narg, char **arg)
{
if (narg < 2) utils::missing_cmd_args(FLERR, "fitpod", error);
@ -241,7 +245,7 @@ void FitPOD::get_exyz_files(std::vector<std::string>& files, const std::string &
{
auto allfiles = platform::list_directory(datapath);
std::sort(allfiles.begin(), allfiles.end());
for (auto fname : allfiles) {
for (const auto &fname : allfiles) {
if (utils::strmatch(fname, fmt::format(".*\\.{}$", extension)))
files.push_back(datapath + platform::filepathsep + fname);
}
@ -252,7 +256,6 @@ int FitPOD::get_number_atom_exyz(std::vector<int>& num_atom, int& num_atom_sum,
std::string filename = file;
FILE *fp;
if (comm->me == 0) {
fp = utils::open_potential(filename,lmp,nullptr);
if (fp == nullptr)
error->one(FLERR,"Cannot open POD coefficient file {}: ", filename, utils::getsyserror());
@ -469,7 +472,7 @@ void FitPOD::get_data(datastruct &data, std::vector<std::string> species)
utils::logmesg(lmp, "{}\n {:^{}} | number of configurations | number of atoms\n{}\n",
sepline, "data file", maxname, sepline);
int i = 0;
for (auto fname : data.data_files) {
for (const auto &fname : data.data_files) {
std::string filename = fname.substr(data.data_path.size()+1);
data.filenames.push_back(filename);
if (comm->me == 0)
@ -695,7 +698,7 @@ void FitPOD::select_data(datastruct &newdata, datastruct data)
data.copydatainfo(newdata);
size_t maxname = 9;
for (auto fname : data.data_files) maxname = MAX(maxname,fname.size());
for (const auto &fname : data.data_files) maxname = MAX(maxname,fname.size());
maxname -= data.data_path.size()+1;
if (comm->me == 0)
@ -872,7 +875,7 @@ int FitPOD::podfullneighborlist(double *y, int *alist, int *neighlist, int *numn
return nn;
}
void FitPOD::allocate_memory(datastruct data)
void FitPOD::allocate_memory(const datastruct &data)
{
int nd = podptr->pod.nd;
memory->create(desc.gd, nd, "fitpod:desc_gd");
@ -989,7 +992,7 @@ void FitPOD::allocate_memory(datastruct data)
}
}
void FitPOD::linear_descriptors(datastruct data, int ci)
void FitPOD::linear_descriptors(const datastruct &data, int ci)
{
int dim = 3;
int nd1 = podptr->pod.nd1;
@ -1020,7 +1023,7 @@ void FitPOD::linear_descriptors(datastruct data, int ci)
}
void FitPOD::quadratic_descriptors(datastruct data, int ci)
void FitPOD::quadratic_descriptors(const datastruct &data, int ci)
{
int dim = 3;
int natom = data.num_atom[ci];
@ -1101,7 +1104,7 @@ void FitPOD::quadratic_descriptors(datastruct data, int ci)
desc.gdd[dim*natom*nd1234+i] = desc.gdd[dim*natom*nd1234+i]/(natom);
}
void FitPOD::cubic_descriptors(datastruct data, int ci)
void FitPOD::cubic_descriptors(const datastruct &data, int ci)
{
int dim = 3;
int natom = data.num_atom[ci];
@ -1168,7 +1171,7 @@ void FitPOD::cubic_descriptors(datastruct data, int ci)
desc.gdd[i] = desc.gdd[i]/(natom*natom);
}
void FitPOD::least_squares_matrix(datastruct data, int ci)
void FitPOD::least_squares_matrix(const datastruct &data, int ci)
{
int dim = 3;
int natom = data.num_atom[ci];
@ -1214,7 +1217,7 @@ void FitPOD::least_squares_matrix(datastruct data, int ci)
}
void FitPOD::least_squares_fit(datastruct data)
void FitPOD::least_squares_fit(const datastruct &data)
{
if (comm->me == 0)
utils::logmesg(lmp, "**************** Begin of Least-Squares Fitting ****************\n");
@ -1308,7 +1311,7 @@ void FitPOD::least_squares_fit(datastruct data)
}
}
double FitPOD::energyforce_calculation(double *force, double *coeff, datastruct data, int ci)
double FitPOD::energyforce_calculation(double *force, double *coeff, const datastruct &data, int ci)
{
int dim = 3;
int *pbc = podptr->pod.pbc;
@ -1348,7 +1351,7 @@ double FitPOD::energyforce_calculation(double *force, double *coeff, datastruct
return energy;
}
void FitPOD::print_analysis(datastruct data, double *outarray, double *errors)
void FitPOD::print_analysis(const datastruct &data, double *outarray, double *errors)
{
int nfiles = data.data_files.size(); // number of files
int lm = 10;
@ -1422,7 +1425,7 @@ void FitPOD::print_analysis(datastruct data, double *outarray, double *errors)
fclose(fp_analysis);
}
void FitPOD::error_analysis(datastruct data, double *coeff)
void FitPOD::error_analysis(const datastruct &data, double *coeff)
{
int dim = 3;
double energy;
@ -1537,7 +1540,10 @@ void FitPOD::error_analysis(datastruct data, double *coeff)
nforceall += nforce;
ci += 1;
}
int q = file + 1;
if (nconfigs == 0) nconfigs = 1;
if (nforceall == 0) nforceall = 1;
errors[0 + 4*q] = emae/nconfigs;
errors[1 + 4*q] = sqrt(essr/nconfigs);
errors[2 + 4*q] = fmae/nforceall;
@ -1550,6 +1556,8 @@ void FitPOD::error_analysis(datastruct data, double *coeff)
errors[3] += fssr;
}
if (nc == 0) nc = 1;
if (nf == 0) nf = 1;
errors[0] = errors[0]/nc;
errors[1] = sqrt(errors[1]/nc);
errors[2] = errors[2]/nf;
@ -1561,7 +1569,7 @@ void FitPOD::error_analysis(datastruct data, double *coeff)
}
}
void FitPOD::energyforce_calculation(datastruct data, double *coeff)
void FitPOD::energyforce_calculation(const datastruct &data, double *coeff)
{
int dim = 3;
double energy;
@ -1816,4 +1824,3 @@ void FitPOD::triclinic_lattice_conversion(double *a, double *b, double *c, doubl
b[0] = bx; b[1] = by; b[2] = 0.0;
c[0] = cx; c[1] = cy; c[2] = cz;
}

View File

@ -26,9 +26,11 @@ CommandStyle(fitpod,FitPOD);
namespace LAMMPS_NS {
class FitPOD : public Command {
private:
public:
FitPOD(LAMMPS *);
void command(int, char **) override;
private:
struct datastruct {
std::string file_format = "extxyz";
std::string file_extension = "xyz";
@ -120,9 +122,6 @@ public:
neighborstruct nb;
class MLPOD *podptr;
FitPOD(LAMMPS *lmp) : Command(lmp) {}
// functions for collecting/collating arrays
void print_matrix(const char *desc, int m, int n, int *a, int lda);
@ -148,7 +147,7 @@ public:
// functions for reading input files and fitting
void command(int, char **) override;
int read_data_file(double *fitting_weights, std::string &file_format, std::string &file_extension,
std::string &test_path, std::string &training_path, std::string &filenametag, const std::string &data_file);
void get_exyz_files(std::vector<std::string> &, const std::string &, const std::string &);
@ -166,16 +165,16 @@ public:
int podneighborlist(int *neighlist, int *numneigh, double *r, double rcutsq, int nx, int N, int dim);
int podfullneighborlist(double *y, int *alist, int *neighlist, int *numneigh, int *numneighsum,
double *x, double *a1, double *a2, double *a3, double rcut, int *pbc, int nx);
void allocate_memory(datastruct data);
void linear_descriptors(datastruct data, int ci);
void quadratic_descriptors(datastruct data, int ci);
void cubic_descriptors(datastruct data, int ci);
void least_squares_matrix(datastruct data, int ci);
void least_squares_fit(datastruct data);
void print_analysis(datastruct data, double *outarray, double *errors);
void error_analysis(datastruct data, double *coeff);
double energyforce_calculation(double *force, double *coeff, datastruct data, int ci);
void energyforce_calculation(datastruct data, double *coeff);
void allocate_memory(const datastruct &data);
void linear_descriptors(const datastruct &data, int ci);
void quadratic_descriptors(const datastruct &data, int ci);
void cubic_descriptors(const datastruct &data, int ci);
void least_squares_matrix(const datastruct &data, int ci);
void least_squares_fit(const datastruct &data);
void print_analysis(const datastruct &data, double *outarray, double *errors);
void error_analysis(const datastruct &data, double *coeff);
double energyforce_calculation(double *force, double *coeff, const datastruct &data, int ci);
void energyforce_calculation(const datastruct &data, double *coeff);
};
} // namespace LAMMPS_NS

View File

@ -695,8 +695,9 @@ void MLPOD::read_coeff_file(const std::string &coeff_file)
/*********************************************************************************************************/
void MLPOD::linear_descriptors(double *gd, double *efatom, double *y, double *tmpmem, int *atomtype,
int *alist, int *pairlist, int *pairnum, int *pairnumsum, int *tmpint, int natom, int Nij)
void MLPOD::linear_descriptors(double *gd, double *efatom, double *y, double *tmpmem,
int *atomtype, int *alist, int *pairlist, int * /*pairnum*/,
int *pairnumsum, int *tmpint, int natom, int Nij)
{
int dim = 3;
int nelements = pod.nelements;
@ -1641,39 +1642,37 @@ void MLPOD::InitSnap()
double rcutmax = 0.0;
for (int ielem = 0; ielem < ntypes; ielem++)
rcutmax = MAX(2.0*elemradius[ielem]*rcutfac,rcutmax);
rcutmax = MAX(2.0*elemradius[ielem]*rcutfac,rcutmax);
snapSetup(twojmax, ntypes);
//TemplateCopytoDevice(&sna.radelem[1], elemradius, ntypes, backend);
//TemplateCopytoDevice(&sna.wjelem[1], elemweight, ntypes, backend);
for (int i=0; i<ntypes; i++) {
sna.radelem[1+i] = elemradius[i];
sna.wjelem[1+i] = elemweight[i];
}
podArrayFill(&sna.map[1], (int) 0, ntypes);
//double cutsq[100];
for (int i=0; i<ntypes; i++)
for (int i=0; i<ntypes; i++) {
for (int j=0; j<ntypes; j++) {
double cut = (elemradius[i] + elemradius[j])*rcutfac;
sna.rcutsq[j+1 + (i+1)*(ntypes+1)] = cut*cut;
}
//TemplateCopytoDevice(sna.rcutsq, cutsq, (ntypes+1)*(ntypes+1), backend);
}
// bzeroflag is currently always 0
#if 0
if (bzeroflag) {
double www = wself*wself*wself;
//double bzero[100];
for (int j = 0; j <= twojmax; j++)
if (bnormflag)
sna.bzero[j] = www;
sna.bzero[j] = www;
else
sna.bzero[j] = www*(j+1);
//TemplateCopytoDevice(sna.bzero, bzero, twojmax+1, backend);
sna.bzero[j] = www*(j+1);
}
#endif
int nelements = ntypes;
if (!chemflag)
nelements = 1;
nelements = 1;
sna.nelements = nelements;
sna.ndoubles = nelements*nelements; // number of multi-element pairs
@ -1902,8 +1901,8 @@ void MLPOD::snapComputeUlist(double *Sr, double *Si, double *dSr, double *dSi, d
};
void MLPOD::snapZeroUarraytot2(double *Stotr, double *Stoti, double wself, int *idxu_block,
int *type, int *map, int *ai, int wselfall_flag, int chemflag, int idxu_max, int nelements,
int twojmax, int inum)
int *type, int *map, int * /*ai*/, int wselfall_flag, int chemflag,
int idxu_max, int nelements, int twojmax, int inum)
{
int N1 = inum;
int N2 = N1*(twojmax+1);
@ -2554,8 +2553,9 @@ void MLPOD::pod1body(double *eatom, int *atomtype, int nelements, int natom)
eatom[i + natom*(m-1)] = (atomtype[i] == m) ? 1.0 : 0.0;
}
void MLPOD::pod3body(double *eatom, double *yij, double *e2ij, double *tmpmem, int *elemindex, int *pairnumsum,
int *idxi, int *ti, int *tj, int nrbf, int nabf, int nelements, int natom, int Nij)
void MLPOD::pod3body(double *eatom, double *yij, double *e2ij, double *tmpmem, int *elemindex,
int *pairnumsum, int * /*idxi*/, int *ti, int *tj, int nrbf, int nabf,
int nelements, int natom, int Nij)
{
int dim = 3, nabf1 = nabf + 1;
int nelements2 = nelements*(nelements+1)/2;
@ -2599,7 +2599,6 @@ void MLPOD::pod3body(double *eatom, double *yij, double *e2ij, double *tmpmem, i
costhe = xdot/(rij*rik);
costhe = costhe > 1.0 ? 1.0 : costhe;
costhe = costhe < -1.0 ? -1.0 : costhe;
xdot = costhe*(rij*rik);
theta = acos(costhe);
for (int p=0; p <nabf1; p++)
@ -3070,8 +3069,8 @@ double MLPOD::calculate_energy(double *energycoeff, double *forcecoeff, double *
return energy;
}
void MLPOD::pod2body_force(double *force, double *fij, double *coeff2, int *ai, int *aj,
int *ti, int *tj, int *elemindex, int nelements, int nbf, int natom, int N)
void MLPOD::pod2body_force(double *force, double *fij, double *coeff2, int *ai, int *aj, int *ti,
int *tj, int *elemindex, int nelements, int nbf, int /*natom*/, int N)
{
int nelements2 = nelements*(nelements+1)/2;
for (int n=0; n<N; n++) {
@ -3236,7 +3235,7 @@ void MLPOD::pod3body_force(double *force, double *yij, double *e2ij, double *f2i
}
void MLPOD::snapTallyForce(double *force, double *dbdr, double *coeff4,
int *ai, int *aj, int *ti, int ijnum, int ncoeff, int ntype)
int *ai, int *aj, int *ti, int ijnum, int ncoeff, int /*ntype*/)
{
int N2 = ijnum*ncoeff;
for (int idx=0; idx<N2; idx++) {
@ -3396,8 +3395,8 @@ double MLPOD::energyforce_calculation(double *force, double *podcoeff, double *e
}
void MLPOD::pod2body_force(double **force, double *fij, double *coeff2, int *ai, int *aj,
int *ti, int *tj, int *elemindex, int nelements, int nbf, int natom, int N)
void MLPOD::pod2body_force(double **force, double *fij, double *coeff2, int *ai, int *aj, int *ti,
int *tj, int *elemindex, int nelements, int nbf, int /*natom*/, int N)
{
int nelements2 = nelements*(nelements+1)/2;
for (int n=0; n<N; n++) {
@ -3559,7 +3558,7 @@ void MLPOD::pod3body_force(double **force, double *yij, double *e2ij, double *f2
}
void MLPOD::snapTallyForce(double **force, double *dbdr, double *coeff4,
int *ai, int *aj, int *ti, int ijnum, int ncoeff, int ntype)
int *ai, int *aj, int *ti, int ijnum, int ncoeff, int /*ntype*/)
{
int N2 = ijnum*ncoeff;
for (int idx=0; idx<N2; idx++) {