tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update kokkospod

Browse Source
This commit is contained in:
exapde
2024-05-18 08:17:48 -04:00
parent feae228329
commit b5fe3d5b06
16 changed files with 1326 additions and 1326 deletions
Download Patch File Download Diff File Expand all files Collapse all files

970
src/KOKKOS/pair_pod_kokkos.cpp
View File

File diff suppressed because it is too large Load Diff

140
src/KOKKOS/pair_pod_kokkos.h
View File

@ -32,7 +32,7 @@ namespace LAMMPS_NS {
template<class DeviceType>
class PairPODKokkos : public PairPOD {
public:
public:
typedef DeviceType device_type;
typedef ArrayTypes<DeviceType> AT;
@ -43,7 +43,7 @@ class PairPODKokkos : public PairPOD {
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
//protected:
int inum, maxneigh;
int host_flag;
@ -53,7 +53,7 @@ class PairPODKokkos : public PairPOD {
typename AT::t_neighbors_2d d_neighbors;
typename AT::t_int_1d d_ilist;
typename AT::t_int_1d d_numneigh;
typename AT::t_int_1d d_numneigh;
// typename AT::t_int_1d_randomread d_ilist;
// typename AT::t_int_1d_randomread d_numneigh;
@ -74,34 +74,34 @@ class PairPODKokkos : public PairPOD {
friend void pair_virial_fdotr_compute<PairPODKokkos>(PairPODKokkos*);
void grow(int, int);
void grow(int, int);
void copy_from_pod_class(EAPOD *podptr);
void divideInterval(int *intervals, int N, int M);
int calculateNumberOfIntervals(int N, int intervalSize);
int calculateNumberOfIntervals(int N, int intervalSize);
void grow_atoms(int Ni);
void grow_pairs(int Nij);
void allocate() override;
double memory_usage() override;
typedef Kokkos::View<int*, DeviceType> t_pod_1i;
typedef Kokkos::View<int*, DeviceType> t_pod_1i;
typedef Kokkos::View<double*, DeviceType> t_pod_1d;
// typedef Kokkos::View<int**, DeviceType> t_pod_2i;
// typedef Kokkos::View<double**, DeviceType> t_pod_2d;
// typedef Kokkos::View<double**[3], DeviceType> t_pod_3d3;
int atomBlockSize; // size of each atom block
int nAtomBlocks; // number of atoms blocks
int atomBlocks[101]; // atom blocks
double comptime[100];
int timing;
int ni; // number of atoms i in the current atom block
int nij; // number of pairs (i,j) in the current atom block
int ni; // number of atoms i in the current atom block
int nij; // number of pairs (i,j) in the current atom block
int nimax; // maximum number of atoms i
int nijmax; // maximum number of pairs (i,j)
int nelements; // number of elements
int nijmax; // maximum number of pairs (i,j)
int nelements; // number of elements
int onebody; // one-body descriptors
int besseldegree; // degree of Bessel functions
int inversedegree; // degree of inverse functions
@ -110,116 +110,116 @@ class PairPODKokkos : public PairPOD {
int ns; // number of snapshots for radial basis functions
int nl1, nl2, nl3, nl4, nl23, nl33, nl34, nl44, nl; // number of local descriptors
int nrbf2, nrbf3, nrbf4, nrbfmax; // number of radial basis functions
int nabf3, nabf4; // number of angular basis functions
int nabf3, nabf4; // number of angular basis functions
int K3, K4, Q4; // number of monomials
// environmental variables
int nClusters; // number of environment clusters
int nComponents; // number of principal components
int Mdesc; // number of base descriptors
int Mdesc; // number of base descriptors
double rin; // inner cut-off radius
double rcut; // outer cut-off radius
double rmax; // rcut - rin
double rmax; // rcut - rin
double rcutsq;
t_pod_1d rij; // (xj - xi) for all pairs (I, J)
t_pod_1d fij; // force for all pairs (I, J)
t_pod_1d ei; // energy for each atom I
t_pod_1i typeai; // types of atoms I only
t_pod_1i numij; // number of pairs (I, J) for each atom I
t_pod_1i numij; // number of pairs (I, J) for each atom I
t_pod_1i idxi; // storing linear indices of atom I for all pairs (I, J)
t_pod_1i ai; // IDs of atoms I for all pairs (I, J)
t_pod_1i aj; // IDs of atoms J for all pairs (I, J)
t_pod_1i ti; // types of atoms I for all pairs (I, J)
t_pod_1i tj; // types of atoms J for all pairs (I, J)
t_pod_1i tj; // types of atoms J for all pairs (I, J)
t_pod_1d besselparams;
t_pod_1d Phi; // eigenvectors matrix ns x ns
t_pod_1d rbf; // radial basis functions nij x nrbfmax
t_pod_1d rbfx; // x-derivatives of radial basis functions nij x nrbfmax
t_pod_1d rbf; // radial basis functions nij x nrbfmax
t_pod_1d rbfx; // x-derivatives of radial basis functions nij x nrbfmax
t_pod_1d rbfy; // y-derivatives of radial basis functions nij x nrbfmax
t_pod_1d rbfz; // z-derivatives of radial basis functions nij x nrbfmax
t_pod_1d rbfz; // z-derivatives of radial basis functions nij x nrbfmax
t_pod_1d abf; // angular basis functions nij x K3
t_pod_1d abfx; // x-derivatives of angular basis functions nij x K3
t_pod_1d abfy; // y-derivatives of angular basis functions nij x K3
t_pod_1d abfy; // y-derivatives of angular basis functions nij x K3
t_pod_1d abfz; // z-derivatives of angular basis functions nij x K3
t_pod_1d sumU; // sum of radial basis functions ni x K3 x nrbfmax x nelements
t_pod_1d forcecoeff; // force coefficients ni x K3 x nrbfmax x nelements
t_pod_1d Proj; // PCA Projection matrix
t_pod_1d Centroids; // centroids of the clusters
t_pod_1d Centroids; // centroids of the clusters
t_pod_1d bd; // base descriptors ni x Mdesc
t_pod_1d cb; // force coefficients for base descriptors ni x Mdesc
t_pod_1d pd; // environment probability descriptors ni x (1 + nComponents + 3*nClusters)
t_pod_1d coefficients; // coefficients nCoeffPerElement x nelements
t_pod_1i pq3, pn3, pc3; // arrays to compute 3-body angular basis functions
t_pod_1i pa4, pb4, pc4; // arrays to compute 4-body angular basis functions
t_pod_1i pa4, pb4, pc4; // arrays to compute 4-body angular basis functions
t_pod_1i ind33l, ind33r; // nl33
t_pod_1i ind34l, ind34r; // nl34
t_pod_1i ind44l, ind44r; // nl44
t_pod_1i elemindex;
t_pod_1i elemindex;
void set_array_to_zero(t_pod_1d a, int N);
int NeighborCount(t_pod_1i, double, int, int);
void NeighborList(t_pod_1d l_rij, t_pod_1i l_numij, t_pod_1i l_typeai, t_pod_1i l_idxi,
void NeighborList(t_pod_1d l_rij, t_pod_1i l_numij, t_pod_1i l_typeai, t_pod_1i l_idxi,
t_pod_1i l_ai, t_pod_1i l_aj, t_pod_1i l_ti, t_pod_1i l_tj, double l_rcutsq, int gi1, int Ni);
void radialbasis(t_pod_1d rbft, t_pod_1d rbftx, t_pod_1d rbfty, t_pod_1d rbftz,
t_pod_1d rij, t_pod_1d l_besselparams, double l_rin, double l_rmax, int l_besseldegree,
int l_inversedegree, int l_nbesselpars, int l_ns, int Nij);
void matrixMultiply(t_pod_1d a, t_pod_1d b, t_pod_1d c, int r1, int c1, int c2);
void radialbasis(t_pod_1d rbft, t_pod_1d rbftx, t_pod_1d rbfty, t_pod_1d rbftz,
t_pod_1d rij, t_pod_1d l_besselparams, double l_rin, double l_rmax, int l_besseldegree,
int l_inversedegree, int l_nbesselpars, int l_ns, int Nij);
void matrixMultiply(t_pod_1d a, t_pod_1d b, t_pod_1d c, int r1, int c1, int c2);
void angularbasis(t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz,
t_pod_1d l_rij, t_pod_1i l_pq3, int l_K3, int N);
void radialangularsum(t_pod_1d l_sumU, t_pod_1d l_rbf, t_pod_1d l_abf, t_pod_1i l_tj,
t_pod_1d l_rij, t_pod_1i l_pq3, int l_K3, int N);
void radialangularsum(t_pod_1d l_sumU, t_pod_1d l_rbf, t_pod_1d l_abf, t_pod_1i l_tj,
t_pod_1i l_numij, const int l_nelements, const int l_nrbf3, const int l_K3, const int Ni, const int Nij);
void twobodydesc(t_pod_1d d2, t_pod_1d l_rbf, t_pod_1i l_idxi, t_pod_1i l_tj, int l_nrbf2, const int Ni, const int Nij);
void threebodydesc(t_pod_1d d3, t_pod_1d l_sumU, t_pod_1i l_pc3, t_pod_1i l_pn3,
void twobodydesc(t_pod_1d d2, t_pod_1d l_rbf, t_pod_1i l_idxi, t_pod_1i l_tj, int l_nrbf2, const int Ni, const int Nij);
void threebodydesc(t_pod_1d d3, t_pod_1d l_sumU, t_pod_1i l_pc3, t_pod_1i l_pn3,
int l_nelements, int l_nrbf3, int l_nabf3, int l_K3, const int Ni);
void fourbodydesc(t_pod_1d d4, t_pod_1d l_sumU, t_pod_1i l_pa4, t_pod_1i l_pb4, t_pod_1i l_pc4,
void fourbodydesc(t_pod_1d d4, t_pod_1d l_sumU, t_pod_1i l_pa4, t_pod_1i l_pb4, t_pod_1i l_pc4,
int l_nelements, int l_nrbf3, int l_nrbf4, int l_nabf4, int l_K3, int l_Q4, int Ni);
void crossdesc(t_pod_1d d12, t_pod_1d d1, t_pod_1d d2, t_pod_1i ind1, t_pod_1i ind2, int n12, int Ni);
void crossdesc_reduction(t_pod_1d cb1, t_pod_1d cb2, t_pod_1d c12, t_pod_1d d1,
t_pod_1d d2, t_pod_1i ind1, t_pod_1i ind2, int n12, int Ni);
void crossdesc_reduction(t_pod_1d cb1, t_pod_1d cb2, t_pod_1d c12, t_pod_1d d1,
t_pod_1d d2, t_pod_1i ind1, t_pod_1i ind2, int n12, int Ni);
void blockatom_base_descriptors(t_pod_1d bd, int Ni, int Nij);
void blockatom_base_coefficients(t_pod_1d ei, t_pod_1d cb, t_pod_1d B, int Ni);
void blockatom_environment_descriptors(t_pod_1d ei, t_pod_1d cb, t_pod_1d B, int Ni);
void twobody_forces(t_pod_1d fij, t_pod_1d cb2, t_pod_1d l_rbfx, t_pod_1d l_rbfy, t_pod_1d l_rbfz,
void blockatom_environment_descriptors(t_pod_1d ei, t_pod_1d cb, t_pod_1d B, int Ni);
void twobody_forces(t_pod_1d fij, t_pod_1d cb2, t_pod_1d l_rbfx, t_pod_1d l_rbfy, t_pod_1d l_rbfz,
t_pod_1i l_idxi, t_pod_1i l_tj, int l_nrbf2, const int Ni, const int Nij);
void threebody_forces(t_pod_1d fij, t_pod_1d cb3, t_pod_1d l_rbf, t_pod_1d l_rbfx,
t_pod_1d l_rbfy, t_pod_1d l_rbfz, t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz,
t_pod_1d l_sumU, t_pod_1i l_idxi, t_pod_1i l_tj, t_pod_1i l_pc3, t_pod_1i l_pn3, t_pod_1i l_elemindex,
void threebody_forces(t_pod_1d fij, t_pod_1d cb3, t_pod_1d l_rbf, t_pod_1d l_rbfx,
t_pod_1d l_rbfy, t_pod_1d l_rbfz, t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz,
t_pod_1d l_sumU, t_pod_1i l_idxi, t_pod_1i l_tj, t_pod_1i l_pc3, t_pod_1i l_pn3, t_pod_1i l_elemindex,
int l_nelements, int l_nrbf3, int l_nabf3, int l_K3, int Ni, int Nij);
void fourbody_forces(t_pod_1d fij, t_pod_1d cb4, t_pod_1d l_rbf, t_pod_1d l_rbfx, t_pod_1d l_rbfy,
t_pod_1d l_rbfz, t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz, t_pod_1d l_sumU,
t_pod_1i l_idxi, t_pod_1i l_tj, t_pod_1i l_pa4, t_pod_1i l_pb4, t_pod_1i l_pc4, t_pod_1i l_elemindex,
void fourbody_forces(t_pod_1d fij, t_pod_1d cb4, t_pod_1d l_rbf, t_pod_1d l_rbfx, t_pod_1d l_rbfy,
t_pod_1d l_rbfz, t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz, t_pod_1d l_sumU,
t_pod_1i l_idxi, t_pod_1i l_tj, t_pod_1i l_pa4, t_pod_1i l_pb4, t_pod_1i l_pc4, t_pod_1i l_elemindex,
int l_nelements, int l_nrbf3, int l_nrbf4, int l_nabf4, int l_K3, int l_Q4, int Ni, int Nij);
void threebody_forcecoeff(t_pod_1d fb3, t_pod_1d cb3, t_pod_1d l_sumU, t_pod_1i l_pc3,
void threebody_forcecoeff(t_pod_1d fb3, t_pod_1d cb3, t_pod_1d l_sumU, t_pod_1i l_pc3,
t_pod_1i l_pn3, t_pod_1i l_elemindex, int l_nelements, int l_nrbf3, int l_nabf3, int l_K3, int Ni);
void fourbody_forcecoeff(t_pod_1d fb4, t_pod_1d cb4, t_pod_1d l_sumU, t_pod_1i l_pa4,
void fourbody_forcecoeff(t_pod_1d fb4, t_pod_1d cb4, t_pod_1d l_sumU, t_pod_1i l_pa4,
t_pod_1i l_pb4, t_pod_1i l_pc4, int l_nelements, int l_nrbf3, int l_nrbf4, int l_nabf4, int l_K3, int l_Q4, int Ni);
void allbody_forces(t_pod_1d fij, t_pod_1d l_forcecoeff, t_pod_1d l_rbf, t_pod_1d l_rbfx,
t_pod_1d l_rbfy, t_pod_1d l_rbfz, t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz,
void allbody_forces(t_pod_1d fij, t_pod_1d l_forcecoeff, t_pod_1d l_rbf, t_pod_1d l_rbfx,
t_pod_1d l_rbfy, t_pod_1d l_rbfz, t_pod_1d l_abf, t_pod_1d l_abfx, t_pod_1d l_abfy, t_pod_1d l_abfz,
t_pod_1i l_idxi, t_pod_1i l_tj, int l_nelements, int l_nrbf3, int l_nabf3, int l_K3, int Nij);
void blockatom_energyforce(t_pod_1d l_ei, t_pod_1d l_fij, int Ni, int Nij);
void tallyenergy(t_pod_1d l_ei, int istart, int Ni);
void tallyforce(t_pod_1d l_fij, t_pod_1i l_ai, t_pod_1i l_aj, int Nij);
void tallystress(t_pod_1d l_fij, t_pod_1d l_rij, t_pod_1i l_ai, t_pod_1i l_aj, int Nij);
void tallyforce(t_pod_1d l_fij, t_pod_1i l_ai, t_pod_1i l_aj, int Nij);
void tallystress(t_pod_1d l_fij, t_pod_1d l_rij, t_pod_1i l_ai, t_pod_1i l_aj, int Nij);
void savematrix2binfile(std::string filename, t_pod_1d d_A, int nrows, int ncols);
void saveintmatrix2binfile(std::string filename, t_pod_1i d_A, int nrows, int ncols);
void savedatafordebugging();

62
src/ML-POD/compute_pod_atom.cpp
View File

@ -34,33 +34,33 @@ enum{SCALAR,VECTOR,ARRAY};
ComputePODAtom::ComputePODAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), list(nullptr), map(nullptr), pod(nullptr), elements(nullptr)
{
{
int nargmin = 7;
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
if (comm->nprocs > 1) error->all(FLERR, "compute command does not support multi processors");
std::string pod_file = std::string(arg[3]); // pod input file
std::string coeff_file = ""; // coefficient input file
std::string proj_file = std::string(arg[4]); // coefficient input file
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
int ntypes = atom->ntypes;
memory->create(map, ntypes + 1, "compute_pod_global:map");
map_element2type(narg - 6, arg + 6, podptr->nelements);
//size_array_rows = 1 + 3*atom->natoms;
map_element2type(narg - 6, arg + 6, podptr->nelements);
//size_array_rows = 1 + 3*atom->natoms;
//size_array_cols = podptr->nCoeffAll;
cutmax = podptr->rcut;
nmax = 0;
nmax = 0;
nijmax = 0;
pod = nullptr;
elements = nullptr;
elements = nullptr;
size_peratom_cols = podptr->Mdesc * podptr->nClusters;
peratom_flag = 1;
}
@ -120,23 +120,23 @@ void ComputePODAtom::compute_peratom()
for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) {
pod[i][icoeff] = 0.0;
}
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list);
double **x = atom->x;
int **firstneigh = list->firstneigh;
int *numneigh = list->numneigh;
int *type = atom->type;
int *ilist = list->ilist;
int inum = list->inum;
int inum = list->inum;
int nClusters = podptr->nClusters;
int Mdesc = podptr->Mdesc;
double rcutsq = podptr->rcut*podptr->rcut;
for (int ii = 0; ii < inum; ii++) {
int i = ilist[ii];
int i = ilist[ii];
int jnum = numneigh[i];
// allocate temporary memory
@ -145,32 +145,32 @@ void ComputePODAtom::compute_peratom()
podptr->free_temp_memory();
podptr->allocate_temp_memory(nijmax);
}
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
ti = &podptr->tmpint[2*nijmax];
tj = &podptr->tmpint[3*nijmax];
// get neighbor list for atom i
lammpsNeighborList(x, firstneigh, atom->tag, type, numneigh, rcutsq, i);
if (nij > 0) {
// peratom base descriptors
double *bd = &podptr->bd[0];
double *bdd = &podptr->bdd[0];
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
if (nClusters>1) {
// peratom env descriptors
double *pd = &podptr->pd[0];
double *pdd = &podptr->pdd[0];
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
for (int k = 0; k < nClusters; k++)
for (int m = 0; m < Mdesc; m++) {
int mk = m + Mdesc*k;
pod[i][mk] = pd[k]*bd[m];
pod[i][mk] = pd[k]*bd[m];
// for (int n=0; n<nij; n++) {
// int ain = 3*ai[n];
// int ajn = 3*aj[n];
@ -182,7 +182,7 @@ void ComputePODAtom::compute_peratom()
// pod[1 + ajn][imk] -= bdd[0 + nm]*pd[k] + bd[m]*pdd[0+nk];
// pod[2 + ajn][imk] -= bdd[1 + nm]*pd[k] + bd[m]*pdd[1+nk];
// pod[3 + ajn][imk] -= bdd[2 + nm]*pd[k] + bd[m]*pdd[2+nk];
// }
// }
}
}
else {
@ -198,11 +198,11 @@ void ComputePODAtom::compute_peratom()
// pod[1 + ajn][im] -= bdd[0 + nm];
// pod[2 + ajn][im] -= bdd[1 + nm];
// pod[3 + ajn][im] -= bdd[2 + nm];
// }
// }
}
}
}
}
}
}
}
/* ----------------------------------------------------------------------
@ -217,7 +217,7 @@ double ComputePODAtom::memory_usage()
}
void ComputePODAtom::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
void ComputePODAtom::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
int *numneigh, double rcutsq, int gi)
{
nij = 0;

8
src/ML-POD/compute_pod_atom.h
View File

@ -35,7 +35,7 @@ class ComputePODAtom : public Compute {
void lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtype, int *numneigh,
double rcutsq, int i);
void map_element2type(int narg, char **arg, int nelements);
private:
class NeighList *list;
class EAPOD *podptr;
@ -46,16 +46,16 @@ class ComputePODAtom : public Compute {
int nijmax;
double *tmpmem; // temporary memory
double *rij; // (xj - xi) for all pairs (I, J)
double *rij; // (xj - xi) for all pairs (I, J)
char **elements;
int *map;
int *ai; // IDs of atoms I for all pairs (I, J)
int *aj; // IDs of atoms J for all pairs (I, J)
int *ti; // types of atoms I for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
};
} // namespace LAMMPS_NS
#endif
#endif
#endif

72
src/ML-POD/compute_pod_global.cpp
View File

@ -34,32 +34,32 @@ enum{SCALAR,VECTOR,ARRAY};
ComputePODGlobal::ComputePODGlobal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), list(nullptr), map(nullptr), pod(nullptr), elements(nullptr)
{
{
array_flag = 1;
extarray = 0;
int nargmin = 7;
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
if (comm->nprocs > 1) error->all(FLERR, "compute command does not support multi processors");
std::string pod_file = std::string(arg[3]); // pod input file
std::string coeff_file = ""; // coefficient input file
std::string proj_file = std::string(arg[4]); // coefficient input file
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
int ntypes = atom->ntypes;
memory->create(map, ntypes + 1, "compute_pod_global:map");
map_element2type(narg - 6, arg + 6, podptr->nelements);
size_array_rows = 1 + 3*atom->natoms;
map_element2type(narg - 6, arg + 6, podptr->nelements);
size_array_rows = 1 + 3*atom->natoms;
size_array_cols = podptr->nCoeffAll;
cutmax = podptr->rcut;
nijmax = 0;
pod = nullptr;
elements = nullptr;
elements = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -87,7 +87,7 @@ void ComputePODGlobal::init()
if (modify->get_compute_by_style("pod").size() > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute pod");
// allocate memory for global array
memory->create(pod,size_array_rows,size_array_cols,
"compute_pod_global:pod");
@ -105,7 +105,7 @@ void ComputePODGlobal::init_list(int /*id*/, NeighList *ptr)
void ComputePODGlobal::compute_array()
{
// int ntotal = atom->nlocal + atom->nghost;
// int ntotal = atom->nlocal + atom->nghost;
invoked_peratom = update->ntimestep;
// clear global array
@ -113,23 +113,23 @@ void ComputePODGlobal::compute_array()
for (int irow = 0; irow < size_array_rows; irow++)
for (int icoeff = 0; icoeff < size_array_cols; icoeff++)
pod[irow][icoeff] = 0.0;
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list);
double **x = atom->x;
int **firstneigh = list->firstneigh;
int *numneigh = list->numneigh;
int *type = atom->type;
int *ilist = list->ilist;
int inum = list->inum;
int inum = list->inum;
int nClusters = podptr->nClusters;
int Mdesc = podptr->Mdesc;
int nCoeffPerElement = podptr->nCoeffPerElement;
double rcutsq = podptr->rcut*podptr->rcut;
for (int ii = 0; ii < inum; ii++) {
int i = ilist[ii];
int jnum = numneigh[i];
@ -140,22 +140,22 @@ void ComputePODGlobal::compute_array()
podptr->free_temp_memory();
podptr->allocate_temp_memory(nijmax);
}
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
ti = &podptr->tmpint[2*nijmax];
tj = &podptr->tmpint[3*nijmax];
// get neighbor list for atom i
lammpsNeighborList(x, firstneigh, atom->tag, type, numneigh, rcutsq, i);
if (nij > 0) {
// peratom base descriptors
double *bd = &podptr->bd[0];
double *bdd = &podptr->bdd[0];
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
double *bdd = &podptr->bdd[0];
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
pod[0][nCoeffPerElement*(ti[0]-1)] += 1.0; // one-body descriptor
@ -163,7 +163,7 @@ void ComputePODGlobal::compute_array()
// peratom env descriptors
double *pd = &podptr->pd[0];
double *pdd = &podptr->pdd[0];
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
for (int j = 0; j < nClusters; j++) {
for (int m=0; m<Mdesc; m++) {
@ -171,7 +171,7 @@ void ComputePODGlobal::compute_array()
pod[0][k] += pd[j]*bd[m];
for (int n=0; n<nij; n++) {
int ain = 3*ai[n];
int ajn = 3*aj[n];
int ajn = 3*aj[n];
int nm = 3*n + 3*nij*m;
int nj = 3*n + 3*nij*j;
pod[1+ain][k] += bdd[0 + nm]*pd[j] + bd[m]*pdd[0 + nj];
@ -182,7 +182,7 @@ void ComputePODGlobal::compute_array()
pod[3+ajn][k] -= bdd[2 + nm]*pd[j] + bd[m]*pdd[2 + nj];
}
}
}
}
}
else {
for (int m=0; m<Mdesc; m++) {
@ -190,7 +190,7 @@ void ComputePODGlobal::compute_array()
pod[0][k] += bd[m];
for (int n=0; n<nij; n++) {
int ain = 3*ai[n];
int ajn = 3*aj[n];
int ajn = 3*aj[n];
int nm = 3*n + 3*nij*m;
pod[1+ain][k] += bdd[0 + nm];
pod[2+ain][k] += bdd[1 + nm];
@ -199,10 +199,10 @@ void ComputePODGlobal::compute_array()
pod[2+ajn][k] -= bdd[1 + nm];
pod[3+ajn][k] -= bdd[2 + nm];
}
}
}
}
}
}
}
}
}
/* ----------------------------------------------------------------------
@ -217,7 +217,7 @@ double ComputePODGlobal::memory_usage()
}
void ComputePODGlobal::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
void ComputePODGlobal::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
int *numneigh, double rcutsq, int gi)
{
nij = 0;

6
src/ML-POD/compute_pod_global.h
View File

@ -35,7 +35,7 @@ class ComputePODGlobal : public Compute {
void lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtype, int *numneigh,
double rcutsq, int i);
void map_element2type(int narg, char **arg, int nelements);
private:
class NeighList *list;
class EAPOD *podptr;
@ -43,7 +43,7 @@ class ComputePODGlobal : public Compute {
double cutmax;
int nij;
int nijmax;
double *tmpmem; // temporary memory
double *rij; // (xj - xi) for all pairs (I, J)
char **elements;
@ -51,7 +51,7 @@ class ComputePODGlobal : public Compute {
int *ai; // IDs of atoms I for all pairs (I, J)
int *aj; // IDs of atoms J for all pairs (I, J)
int *ti; // types of atoms I for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
};
} // namespace LAMMPS_NS

74
src/ML-POD/compute_pod_local.cpp
View File

@ -34,34 +34,34 @@ enum{SCALAR,VECTOR,ARRAY};
ComputePODLocal::ComputePODLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), list(nullptr), map(nullptr), pod(nullptr), elements(nullptr)
{
{
array_flag = 1;
extarray = 0;
int nargmin = 7;
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
if (comm->nprocs > 1) error->all(FLERR, "compute command does not support multi processors");
std::string pod_file = std::string(arg[3]); // pod input file
std::string coeff_file = ""; // coefficient input file
std::string proj_file = std::string(arg[4]); // coefficient input file
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
int ntypes = atom->ntypes;
memory->create(map, ntypes + 1, "compute_pod_local:map");
map_element2type(narg - 6, arg + 6, podptr->nelements);
map_element2type(narg - 6, arg + 6, podptr->nelements);
int numdesc = podptr->Mdesc * podptr->nClusters;
size_array_rows = 1 + 3*atom->natoms;
size_array_rows = 1 + 3*atom->natoms;
size_array_cols = atom->natoms*numdesc;
cutmax = podptr->rcut;
nijmax = 0;
pod = nullptr;
elements = nullptr;
elements = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -89,7 +89,7 @@ void ComputePODLocal::init()
if (modify->get_compute_by_style("pod").size() > 1 && comm->me == 0)
error->warning(FLERR,"More than one compute pod");
// allocate memory for global array
memory->create(pod,size_array_rows,size_array_cols,
"compute_pod_local:pod");
@ -107,7 +107,7 @@ void ComputePODLocal::init_list(int /*id*/, NeighList *ptr)
void ComputePODLocal::compute_array()
{
// int ntotal = atom->nlocal + atom->nghost;
// int ntotal = atom->nlocal + atom->nghost;
invoked_peratom = update->ntimestep;
// clear global array
@ -115,22 +115,22 @@ void ComputePODLocal::compute_array()
for (int irow = 0; irow < size_array_rows; irow++)
for (int icoeff = 0; icoeff < size_array_cols; icoeff++)
pod[irow][icoeff] = 0.0;
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list);
double **x = atom->x;
int **firstneigh = list->firstneigh;
int *numneigh = list->numneigh;
int *type = atom->type;
int *ilist = list->ilist;
int inum = list->inum;
int inum = list->inum;
int nClusters = podptr->nClusters;
int Mdesc = podptr->Mdesc;
double rcutsq = podptr->rcut*podptr->rcut;
for (int ii = 0; ii < inum; ii++) {
int i = ilist[ii];
int jnum = numneigh[i];
@ -141,32 +141,32 @@ void ComputePODLocal::compute_array()
podptr->free_temp_memory();
podptr->allocate_temp_memory(nijmax);
}
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
ti = &podptr->tmpint[2*nijmax];
tj = &podptr->tmpint[3*nijmax];
// get neighbor list for atom i
lammpsNeighborList(x, firstneigh, atom->tag, type, numneigh, rcutsq, i);
if (nij > 0) {
// peratom base descriptors
double *bd = &podptr->bd[0];
double *bdd = &podptr->bdd[0];
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
double *bdd = &podptr->bdd[0];
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
if (nClusters>1) {
// peratom env descriptors
double *pd = &podptr->pd[0];
double *pdd = &podptr->pdd[0];
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
for (int k = 0; k < nClusters; k++)
for (int m = 0; m < Mdesc; m++) {
int imk = m + Mdesc*k + Mdesc*nClusters*i;
pod[0][imk] = pd[k]*bd[m];
pod[0][imk] = pd[k]*bd[m];
for (int n=0; n<nij; n++) {
int ain = 3*ai[n];
int ajn = 3*aj[n];
@ -178,7 +178,7 @@ void ComputePODLocal::compute_array()
pod[1 + ajn][imk] -= bdd[0 + nm]*pd[k] + bd[m]*pdd[0+nk];
pod[2 + ajn][imk] -= bdd[1 + nm]*pd[k] + bd[m]*pdd[1+nk];
pod[3 + ajn][imk] -= bdd[2 + nm]*pd[k] + bd[m]*pdd[2+nk];
}
}
}
}
else {
@ -195,11 +195,11 @@ void ComputePODLocal::compute_array()
pod[1 + ajn][im] -= bdd[0 + nm];
pod[2 + ajn][im] -= bdd[1 + nm];
pod[3 + ajn][im] -= bdd[2 + nm];
}
}
}
}
}
}
}
}
}
/* ----------------------------------------------------------------------
@ -214,7 +214,7 @@ double ComputePODLocal::memory_usage()
}
void ComputePODLocal::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
void ComputePODLocal::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
int *numneigh, double rcutsq, int gi)
{
nij = 0;
@ -259,7 +259,7 @@ void ComputePODLocal::map_element2type(int narg, char **arg, int nelements)
}
elements = new char*[ntypes];
for (i = 0; i < ntypes; i++) elements[i] = nullptr;
nelements = 0;
map[0] = -1;
for (i = 1; i <= narg; i++) {
@ -275,5 +275,5 @@ void ComputePODLocal::map_element2type(int narg, char **arg, int nelements)
elements[j] = utils::strdup(entry);
nelements++;
}
}
}
}

6
src/ML-POD/compute_pod_local.h
View File

@ -35,7 +35,7 @@ class ComputePODLocal : public Compute {
void lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtype, int *numneigh,
double rcutsq, int i);
void map_element2type(int narg, char **arg, int nelements);
private:
class NeighList *list;
class EAPOD *podptr;
@ -43,7 +43,7 @@ class ComputePODLocal : public Compute {
double cutmax;
int nij;
int nijmax;
double *tmpmem; // temporary memory
double *rij; // (xj - xi) for all pairs (I, J)
char **elements;
@ -51,7 +51,7 @@ class ComputePODLocal : public Compute {
int *ai; // IDs of atoms I for all pairs (I, J)
int *aj; // IDs of atoms J for all pairs (I, J)
int *ti; // types of atoms I for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
};
} // namespace LAMMPS_NS

60
src/ML-POD/compute_podd_atom.cpp
View File

@ -34,30 +34,30 @@ enum{SCALAR,VECTOR,ARRAY};
ComputePODDAtom::ComputePODDAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), list(nullptr), map(nullptr), pod(nullptr), elements(nullptr)
{
{
int nargmin = 7;
if (narg < nargmin) error->all(FLERR, "Illegal compute {} command", style);
if (comm->nprocs > 1) error->all(FLERR, "compute command does not support multi processors");
std::string pod_file = std::string(arg[3]); // pod input file
std::string coeff_file = ""; // coefficient input file
std::string proj_file = std::string(arg[4]); // coefficient input file
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
std::string centroid_file = std::string(arg[5]); // coefficient input file
podptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, centroid_file);
int ntypes = atom->ntypes;
memory->create(map, ntypes + 1, "compute_pod_global:map");
map_element2type(narg - 6, arg + 6, podptr->nelements);
map_element2type(narg - 6, arg + 6, podptr->nelements);
cutmax = podptr->rcut;
nmax = 0;
nmax = 0;
nijmax = 0;
pod = nullptr;
elements = nullptr;
elements = nullptr;
size_peratom_cols = podptr->Mdesc * podptr->nClusters*3*atom->natoms;
peratom_flag = 1;
}
@ -117,23 +117,23 @@ void ComputePODDAtom::compute_peratom()
for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++) {
pod[i][icoeff] = 0.0;
}
// invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list);
double **x = atom->x;
int **firstneigh = list->firstneigh;
int *numneigh = list->numneigh;
int *type = atom->type;
int *ilist = list->ilist;
int inum = list->inum;
int inum = list->inum;
int nClusters = podptr->nClusters;
int Mdesc = podptr->Mdesc;
double rcutsq = podptr->rcut*podptr->rcut;
for (int ii = 0; ii < inum; ii++) {
int i = ilist[ii];
int i = ilist[ii];
int jnum = numneigh[i];
// allocate temporary memory
@ -142,34 +142,34 @@ void ComputePODDAtom::compute_peratom()
podptr->free_temp_memory();
podptr->allocate_temp_memory(nijmax);
}
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
rij = &podptr->tmpmem[0];
tmpmem = &podptr->tmpmem[3*nijmax];
ai = &podptr->tmpint[0];
aj = &podptr->tmpint[nijmax];
ti = &podptr->tmpint[2*nijmax];
tj = &podptr->tmpint[3*nijmax];
// get neighbor list for atom i
lammpsNeighborList(x, firstneigh, atom->tag, type, numneigh, rcutsq, i);
if (nij > 0) {
// peratom base descriptors
double *bd = &podptr->bd[0];
double *bdd = &podptr->bdd[0];
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
podptr->peratombase_descriptors(bd, bdd, rij, tmpmem, ti, tj, nij);
if (nClusters>1) {
// peratom env descriptors
double *pd = &podptr->pd[0];
double *pdd = &podptr->pdd[0];
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
podptr->peratomenvironment_descriptors(pd, pdd, bd, bdd, tmpmem, ti[0] - 1, nij);
for (int n=0; n<nij; n++) {
int ain = 3*ai[n];
int ajn = 3*aj[n];
int ajn = 3*aj[n];
for (int k = 0; k < nClusters; k++) {
for (int m = 0; m < Mdesc; m++) {
int mk = m + Mdesc*k;
int mk = m + Mdesc*k;
int nm = 3*n + 3*nij*m;
int nk = 3*n + 3*nij*k;
pod[i][mk + Mdesc*nClusters*0 + Mdesc*nClusters*ain] += bdd[0 + nm]*pd[k] + bd[m]*pdd[0+nk];
@ -178,7 +178,7 @@ void ComputePODDAtom::compute_peratom()
pod[i][mk + Mdesc*nClusters*0 + Mdesc*nClusters*ajn] -= bdd[0 + nm]*pd[k] + bd[m]*pdd[0+nk];
pod[i][mk + Mdesc*nClusters*1 + Mdesc*nClusters*ajn] -= bdd[1 + nm]*pd[k] + bd[m]*pdd[1+nk];
pod[i][mk + Mdesc*nClusters*2 + Mdesc*nClusters*ajn] -= bdd[2 + nm]*pd[k] + bd[m]*pdd[2+nk];
}
}
}
}
}
@ -194,11 +194,11 @@ void ComputePODDAtom::compute_peratom()
pod[i][m + Mdesc*0 + Mdesc*ajn] -= bdd[0 + nm];
pod[i][m + Mdesc*1 + Mdesc*ajn] -= bdd[1 + nm];
pod[i][m + Mdesc*2 + Mdesc*ajn] -= bdd[2 + nm];
}
}
}
}
}
}
}
}
}
/* ----------------------------------------------------------------------
@ -213,7 +213,7 @@ double ComputePODDAtom::memory_usage()
}
void ComputePODDAtom::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
void ComputePODDAtom::lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtypes,
int *numneigh, double rcutsq, int gi)
{
nij = 0;

8
src/ML-POD/compute_podd_atom.h
View File

@ -35,7 +35,7 @@ class ComputePODDAtom : public Compute {
void lammpsNeighborList(double **x, int **firstneigh, int *atomid, int *atomtype, int *numneigh,
double rcutsq, int i);
void map_element2type(int narg, char **arg, int nelements);
private:
class NeighList *list;
class EAPOD *podptr;
@ -46,16 +46,16 @@ class ComputePODDAtom : public Compute {
int nijmax;
double *tmpmem; // temporary memory
double *rij; // (xj - xi) for all pairs (I, J)
double *rij; // (xj - xi) for all pairs (I, J)
char **elements;
int *map;
int *ai; // IDs of atoms I for all pairs (I, J)
int *aj; // IDs of atoms J for all pairs (I, J)
int *ti; // types of atoms I for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
};
} // namespace LAMMPS_NS
#endif
#endif
#endif

332
src/ML-POD/eapod.cpp
View File

@ -49,7 +49,7 @@ EAPOD::EAPOD(LAMMPS *_lmp, const std::string &pod_file, const std::string &coeff
ind34r = nullptr;
ind44l = nullptr;
ind44r = nullptr;
rin = 0.5;
rcut = 5.0;
nClusters = 1;
@ -101,20 +101,20 @@ EAPOD::EAPOD(LAMMPS *_lmp, const std::string &pod_file, const std::string &coeff
error->all(FLERR,"number of coefficients in the coefficient file is not correct");
}
if (nClusters > 1) {
// read projection matrix file to podstruct
// read projection matrix file to podstruct
if (proj_file != "") {
nproj = read_projection_matrix(proj_file);
if (nproj != nComponents*Mdesc*nelements)
error->all(FLERR,"number of coefficients in the projection file is not correct");
}
// read centroids file to podstruct
if (centroids_file != "") {
ncentroids = read_centroids(centroids_file);
if (ncentroids != nComponents*nClusters*nelements)
error->all(FLERR,"number of coefficients in the projection file is not correct");
}
}
}
}
// destructor
@ -129,7 +129,7 @@ EAPOD::~EAPOD()
memory->destroy(bdd);
memory->destroy(pd);
memory->destroy(pdd);
memory->destroy(coeff);
memory->destroy(coeff);
memory->destroy(tmpmem);
memory->destroy(tmpint);
memory->destroy(pn3);
@ -150,7 +150,7 @@ EAPOD::~EAPOD()
memory->destroy(ind44l);
memory->destroy(ind33r);
memory->destroy(ind34r);
memory->destroy(ind44r);
memory->destroy(ind44r);
}
void EAPOD::read_pod_file(std::string pod_file)
@ -268,7 +268,7 @@ void EAPOD::read_pod_file(std::string pod_file)
if (P4 < P44) error->all(FLERR,"seven-body angular degree must be equal or less than four-body angular degree");
if (P3 > 12) error->all(FLERR,"three-body angular degree must be equal or less than 12");
if (P4 > 6) error->all(FLERR,"four-body angular degree must be equal or less than 6");
if (P4 > 6) error->all(FLERR,"four-body angular degree must be equal or less than 6");
int Ne = nelements;
memory->create(elemindex, Ne*Ne, "elemindex");
@ -416,7 +416,7 @@ void EAPOD::read_pod_file(std::string pod_file)
utils::logmesg(lmp, "number of local descriptors per element for five-body potential: {}\n", nl33);
utils::logmesg(lmp, "number of local descriptors per element for six-body potential: {}\n", nl34);
utils::logmesg(lmp, "number of local descriptors per element for seven-body potential: {}\n", nl44);
utils::logmesg(lmp, "number of local descriptors per element for all potentials: {}\n", nl);
utils::logmesg(lmp, "number of local descriptors per element for all potentials: {}\n", nl);
utils::logmesg(lmp, "number of global descriptors: {}\n", nCoeffAll);
utils::logmesg(lmp, "**************** End of POD Potentials ****************\n\n");
}
@ -702,7 +702,7 @@ void EAPOD::peratombase_descriptors(double *bd1, double *bdd1, double *rij, doub
for (int i=0; i<3*Nj*Mdesc; i++) bdd1[i] = 0.0;
if (Nj == 0) return;
double *d2 = &bd1[0]; // nl2
double *d3 = &bd1[nl2]; // nl3
double *d4 = &bd1[nl2 + nl3]; // nl4
@ -711,8 +711,8 @@ void EAPOD::peratombase_descriptors(double *bd1, double *bdd1, double *rij, doub
double *d34 = &bd1[nl2 + nl3 + nl4 + nl23 + nl33]; // nl34
double *d44 = &bd1[nl2 + nl3 + nl4 + nl23 + nl33 + nl34]; // nl44
double *dd2 = &bdd1[0]; // 3*Nj*nl2
double *dd3 = &bdd1[3*Nj*nl2]; // 3*Nj*nl3
double *dd2 = &bdd1[0]; // 3*Nj*nl2
double *dd3 = &bdd1[3*Nj*nl2]; // 3*Nj*nl3
double *dd4 = &bdd1[3*Nj*(nl2+nl3)]; // 3*Nj*nl4
double *dd23 = &bdd1[3*Nj*(nl2+nl3+nl4)]; // 3*Nj*nl23
double *dd33 = &bdd1[3*Nj*(nl2+nl3+nl4+nl23)]; // 3*Nj*nl33
@ -797,7 +797,7 @@ void EAPOD::peratombase_descriptors(double *bd1, double *bdd1, double *rij, doub
}
}
fourbodydescderiv(d4, dd4, sumU, Ux, Uy, Uz, tj, Nj);
if ((nl34>0) && (Nj>4)) {
crossdesc(d34, d3, d4, ind34l, ind34r, nl34);
crossdescderiv(dd34, d3, d4, dd3, dd4, ind34l, ind34r, nl34, 3*Nj);
@ -812,39 +812,39 @@ void EAPOD::peratombase_descriptors(double *bd1, double *bdd1, double *rij, doub
}
double EAPOD::peratombase_coefficients(double *cb, double *bd, int *ti)
{
{
int nc = nCoeffPerElement*(ti[0]-1);
double ei = coeff[0 + nc];
for (int m=0; m<Mdesc; m++) {
ei += coeff[1 + m + nc]*bd[m];
cb[m] = coeff[1 + m + nc];
}
for (int m=0; m<Mdesc; m++) {
ei += coeff[1 + m + nc]*bd[m];
cb[m] = coeff[1 + m + nc];
}
return ei;
}
double EAPOD::peratom_environment_descriptors(double *cb, double *bd, double *tm, int *ti)
{
double *P = &tm[0]; // nClusters
double *cp = &tm[(nClusters)]; // nClusters
{
double *P = &tm[0]; // nClusters
double *cp = &tm[(nClusters)]; // nClusters
double *D = &tm[(2*nClusters)]; // nClusters
double *pca = &tm[(3*nClusters)]; // nComponents
double *pca = &tm[(3*nClusters)]; // nComponents
double *proj = &Proj[0];
double *cent = &Centroids[0];
double *cent = &Centroids[0];
int typei = ti[0]-1;
for (int k=0; k<nComponents; k++) {
double sum = 0.0;
double sum = 0.0;
for (int m = 0; m < Mdesc; m++) {
sum += proj[k + nComponents*m + nComponents*Mdesc*typei] * bd[m];
}
pca[k] = sum;
pca[k] = sum;
}
for (int j=0; j<nClusters; j++) {
double sum = 1e-20;
double sum = 1e-20;
for (int k = 0; k < nComponents; k++) {
double c = cent[k + j * nComponents + nClusters*nComponents*typei];
double p = pca[k];
@ -852,35 +852,35 @@ double EAPOD::peratom_environment_descriptors(double *cb, double *bd, double *tm
}
D[j] = 1.0 / sum;
}
double sum = 0;
for (int j = 0; j < nClusters; j++) sum += D[j];
double sum = 0;
for (int j = 0; j < nClusters; j++) sum += D[j];
double sumD = sum;
for (int j = 0; j < nClusters; j++) P[j] = D[j]/sum;
for (int j = 0; j < nClusters; j++) P[j] = D[j]/sum;
int nc = nCoeffPerElement*(ti[0]-1);
double ei = coeff[0 + nc];
for (int k = 0; k<nClusters; k++)
for (int m=0; m<Mdesc; m++)
for (int m=0; m<Mdesc; m++)
ei += coeff[1 + m + Mdesc*k + nc]*bd[m]*P[k];
for (int k=0; k<nClusters; k++) {
for (int k=0; k<nClusters; k++) {
double sum = 0;
for (int m = 0; m<Mdesc; m++)
for (int m = 0; m<Mdesc; m++)
sum += coeff[1 + m + k*Mdesc + nc]*bd[m];
cp[k] = sum;
}
}
for (int m = 0; m<Mdesc; m++) {
double sum = 0.0;
for (int k = 0; k<nClusters; k++)
sum += coeff[1 + m + k*Mdesc + nc]*P[k];
double sum = 0.0;
for (int k = 0; k<nClusters; k++)
sum += coeff[1 + m + k*Mdesc + nc]*P[k];
cb[m] = sum;
}
}
for (int m = 0; m<Mdesc; m++) {
double S1 = 1/sumD;
double S2 = sumD*sumD;
double S2 = sumD*sumD;
double sum = 0.0;
for (int j=0; j<nClusters; j++) {
double dP_dB = 0.0;
@ -890,16 +890,16 @@ double EAPOD::peratom_environment_descriptors(double *cb, double *bd, double *tm
double dD_dB = 0.0;
double D2 = 2 * D[k] * D[k];
for (int n = 0; n < nComponents; n++) {
double dD_dpca = D2 * (cent[n + k * nComponents + nClusters*nComponents*typei] - pca[n]);
double dD_dpca = D2 * (cent[n + k * nComponents + nClusters*nComponents*typei] - pca[n]);
dD_dB += dD_dpca * proj[n + m * nComponents + nComponents*Mdesc*typei];
}
}
dP_dB += dP_dD * dD_dB;
}
sum += cp[j]*dP_dB;
}
sum += cp[j]*dP_dB;
}
cb[m] += sum;
}
}
return ei;
}
@ -912,7 +912,7 @@ void EAPOD::twobody_forces(double *fij, double *cb2, double *rbfx, double *rbfy,
int m = idx % nrbf2; // Recalculate n
int i2 = n + Nj * m; // Index of the radial basis function for atom n and RBF m
int i1 = 3*n;
int i1 = 3*n;
double c = cb2[m + nrbf2*(tj[n] - 1)];
fij[0 + i1] += c*rbfx[i2]; // Add the derivative with respect to x to the corresponding descriptor derivative
fij[1 + i1] += c*rbfy[i2]; // Add the derivative with respect to y to the corresponding descriptor derivative
@ -925,20 +925,20 @@ void EAPOD::threebody_forcecoeff(double *fb3, double *cb3, double *sumU)
if (nelements==1) {
for (int m = 0; m < nrbf3; ++m) {
for (int p = 0; p < nabf3; p++) {
double c3 = 2.0 * cb3[p + nabf3*m];
double c3 = 2.0 * cb3[p + nabf3*m];
int n1 = pn3[p];
int n2 = pn3[p + 1];
int nn = n2 - n1;
int idxU = K3 * m;
for (int q = 0; q < nn; q++) {
int k = n1 + q;
fb3[k + idxU] += c3 * pc3[k] * sumU[k + idxU];
int nn = n2 - n1;
int idxU = K3 * m;
for (int q = 0; q < nn; q++) {
int k = n1 + q;
fb3[k + idxU] += c3 * pc3[k] * sumU[k + idxU];
}
}
}
}
else {
int N3 = nabf3 * nrbf3;
int N3 = nabf3 * nrbf3;
for (int m = 0; m < nrbf3; ++m) {
for (int p = 0; p < nabf3; p++) {
int n1 = pn3[p];
@ -947,13 +947,13 @@ void EAPOD::threebody_forcecoeff(double *fb3, double *cb3, double *sumU)
int jmp = p + nabf3*m;
for (int q = 0; q < nn; q++) {
int k = n1 + q; // Combine n1 and q into a single index
int idxU = nelements * k + nelements * K3 * m;
int idxU = nelements * k + nelements * K3 * m;
for (int i1 = 0; i1 < nelements; i1++) {
double tm = pc3[k] * sumU[i1 + idxU];
for (int i2 = i1; i2 < nelements; i2++) {
int em = elemindex[i2 + nelements * i1];
int em = elemindex[i2 + nelements * i1];
double t1 = tm * cb3[jmp + N3*em]; // Ni * nabf3 * nrbf3 * nelements*(nelements+1)/2
fb3[i2 + idxU] += t1; // K3*nrbf3*Ni
fb3[i2 + idxU] += t1; // K3*nrbf3*Ni
fb3[i1 + idxU] += pc3[k] * cb3[jmp + N3*em] * sumU[i2 + idxU];
}
}
@ -967,7 +967,7 @@ void EAPOD::fourbody_forcecoeff(double *fb4, double *cb4, double *sumU)
{
if (nelements==1) {
for (int m = 0; m < nrbf4; ++m) {
int idxU = K3 * m;
int idxU = K3 * m;
for (int p = 0; p < nabf4; p++) {
int n1 = pa4[p];
int n2 = pa4[p + 1];
@ -981,21 +981,21 @@ void EAPOD::fourbody_forcecoeff(double *fb4, double *cb4, double *sumU)
int j3 = idxU + pb4[idxNQ + 2 * Q4];
double c1 = sumU[j1];
double c2 = sumU[j2];
double c3 = sumU[j3];
fb4[j3] += c * c1 * c2;
double c3 = sumU[j3];
fb4[j3] += c * c1 * c2;
fb4[j2] += c * c1 * c3;
fb4[j1] += c * c2 * c3;
fb4[j1] += c * c2 * c3;
}
}
}
}
else {
int N3 = nabf4 * nrbf4;
else {
int N3 = nabf4 * nrbf4;
for (int m = 0; m < nrbf4; ++m) {
for (int p = 0; p < nabf4; p++) {
int n1 = pa4[p];
int n2 = pa4[p + 1];
int nn = n2 - n1;
int nn = n2 - n1;
int jpm = p + nabf4*m;
for (int q = 0; q < nn; q++) {
int c = pc4[n1 + q];
@ -1005,17 +1005,17 @@ void EAPOD::fourbody_forcecoeff(double *fb4, double *cb4, double *sumU)
int idx1 = nelements * j1 + nelements * K3 * m;
int idx2 = nelements * j2 + nelements * K3 * m;
int idx3 = nelements * j3 + nelements * K3 * m;
int k = 0;
for (int i1 = 0; i1 < nelements; i1++) {
int k = 0;
for (int i1 = 0; i1 < nelements; i1++) {
double c1 = sumU[idx1 + i1];
for (int i2 = i1; i2 < nelements; i2++) {
double c2 = sumU[idx2 + i2];
for (int i3 = i2; i3 < nelements; i3++) {
double c3 = sumU[idx3 + i3];
double c2 = sumU[idx2 + i2];
for (int i3 = i2; i3 < nelements; i3++) {
double c3 = sumU[idx3 + i3];
double c4 = c * cb4[jpm + N3*k];
fb4[idx3 + i3] += c4*(c1 * c2);
fb4[idx3 + i3] += c4*(c1 * c2);
fb4[idx2 + i2] += c4*(c1 * c3);
fb4[idx1 + i1] += c4*(c2 * c3);
fb4[idx1 + i1] += c4*(c2 * c3);
k += 1;
}
}
@ -1029,74 +1029,74 @@ void EAPOD::fourbody_forcecoeff(double *fb4, double *cb4, double *sumU)
void EAPOD::allbody_forces(double *fij, double *forcecoeff, double *rbf, double *rbfx, double *rbfy, double *rbfz,
double *abf, double *abfx, double *abfy, double *abfz, int *tj, int Nj)
{
int totalIterations = nrbf3 * Nj;
int totalIterations = nrbf3 * Nj;
for (int idx = 0; idx < totalIterations; ++idx) {
int j = idx / nrbf3; // Derive the original j value
int m = idx % nrbf3; // Derive the original m value
int idxR = m + nrbfmax * j; // Pre-compute the index for rbf
int i2 = tj[j] - 1;
int i2 = tj[j] - 1;
double rbfBase = rbf[idxR];
double rbfxBase = rbfx[idxR];
double rbfyBase = rbfy[idxR];
double rbfzBase = rbfz[idxR];
double fx = 0;
double fy = 0;
double fz = 0;
double fz = 0;
for (int k = 0; k < K3; k++) {
int idxU = nelements * k + nelements * K3 * m;
double fc = forcecoeff[i2 + idxU];
int idxA = k + K3 * j; // Pre-compute the index for abf
double abfA = abf[idxA];
int idxA = k + K3 * j; // Pre-compute the index for abf
double abfA = abf[idxA];
double abfxA = abfx[idxA];
double abfyA = abfy[idxA];
double abfzA = abfz[idxA];
fx += fc * (abfxA * rbfBase + rbfxBase * abfA); // K3*nrbf3*Nij
double abfzA = abfz[idxA];
fx += fc * (abfxA * rbfBase + rbfxBase * abfA); // K3*nrbf3*Nij
fy += fc * (abfyA * rbfBase + rbfyBase * abfA);
fz += fc * (abfzA * rbfBase + rbfzBase * abfA);
fz += fc * (abfzA * rbfBase + rbfzBase * abfA);
}
int baseIdx = 3 * j;
int baseIdx = 3 * j;
fij[baseIdx] += fx;
fij[baseIdx + 1] += fy;
fij[baseIdx + 2] += fz;
}
fij[baseIdx + 2] += fz;
}
}
void EAPOD::allbody_forces(double *fij, double *forcecoeff, double *Ux, double *Uy, double *Uz, int *tj, int Nj)
{
for (int j = 0; j < Nj; ++j) {
int i2 = tj[j] - 1;
int i2 = tj[j] - 1;
double fx = 0;
double fy = 0;
double fz = 0;
for (int m = 0; m < nrbf3; m++)
for (int k = 0; k < K3; k++) {
double fz = 0;
for (int m = 0; m < nrbf3; m++)
for (int k = 0; k < K3; k++) {
double fc = forcecoeff[i2 + nelements * k + nelements * K3 * m];
int idxU = j + Nj*k + Nj * K3 * m; // Pre-compute the index for abf
fx += fc * Ux[idxU]; // K3*nrbf3*Nij
int idxU = j + Nj*k + Nj * K3 * m; // Pre-compute the index for abf
fx += fc * Ux[idxU]; // K3*nrbf3*Nij
fy += fc * Uy[idxU];
fz += fc * Uz[idxU];
fz += fc * Uz[idxU];
}
int baseIdx = 3 * j;
int baseIdx = 3 * j;
fij[baseIdx] += fx;
fij[baseIdx + 1] += fy;
fij[baseIdx + 2] += fz;
}
fij[baseIdx + 2] += fz;
}
}
double EAPOD::peratomenergyforce2(double *fij, double *rij, double *temp,
int *ti, int *tj, int Nj)
{
if (Nj==0) {
return coeff[nCoeffPerElement*(ti[0]-1)];
return coeff[nCoeffPerElement*(ti[0]-1)];
}
int N = 3*Nj;
for (int n=0; n<N; n++) fij[n] = 0.0;
//double *coeff1 = &coeff[nCoeffPerElement*(ti[0]-1)];
double e = 0.0;
double e = 0.0;
for (int i=0; i<Mdesc; i++) bd[i] = 0.0;
double *d2 = &bd[0]; // nl2
double *d3 = &bd[nl2]; // nl3
double *d4 = &bd[nl2 + nl3]; // nl4
@ -1126,7 +1126,7 @@ double EAPOD::peratomenergyforce2(double *fij, double *rij, double *temp,
double *rbfxt = &temp[4*n1 + n5 + 4*n2 + n3]; // Nj*ns
double *rbfyt = &temp[4*n1 + n5 + 4*n2 + 2*n3]; // Nj*ns
double *rbfzt = &temp[4*n1 + n5 + 4*n2 + 3*n3]; // Nj*ns
radialbasis(rbft, rbfxt, rbfyt, rbfzt, rij, besselparams, rin, rcut-rin, pdegree[0], pdegree[1], nbesselpars, Nj);
char chn = 'N';
@ -1144,7 +1144,7 @@ double EAPOD::peratomenergyforce2(double *fij, double *rij, double *temp,
double *abf = &temp[4*n1 + n5 + 4*n2]; // Nj*K3
double *abfx = &temp[4*n1 + n5 + 4*n2 + n4]; // Nj*K3
double *abfy = &temp[4*n1 + n5 + 4*n2 + 2*n4]; // Nj*K3
double *abfz = &temp[4*n1 + n5 + 4*n2 + 3*n4]; // Nj*K3
double *abfz = &temp[4*n1 + n5 + 4*n2 + 3*n4]; // Nj*K3
double *tm = &temp[4*n1 + n5 + 4*n2 + 4*n4]; // 4*K3
angularbasis(abf, abfx, abfy, abfz, rij, tm, pq3, Nj, K3);
@ -1164,7 +1164,7 @@ double EAPOD::peratomenergyforce2(double *fij, double *rij, double *temp,
if ((nl4 > 0) && (Nj>2)) {
fourbodydesc(d4, sumU);
if ((nl34>0) && (Nj>4)) {
crossdesc(d34, d3, d4, ind34l, ind34r, nl34);
}
@ -1174,47 +1174,47 @@ double EAPOD::peratomenergyforce2(double *fij, double *rij, double *temp,
}
}
}
double *cb = &bdd[0];
if (nClusters > 1) {
double *cb = &bdd[0];
if (nClusters > 1) {
e += peratom_environment_descriptors(cb, bd, &temp[4*n1 + n5 + 4*n2], ti);
}
else {
e += peratombase_coefficients(cb, bd, ti);
}
}
double *cb2 = &cb[0]; // nl3
double *cb3 = &cb[nl2]; // nl3
double *cb4 = &cb[(nl2 + nl3)]; // nl4
double *cb33 = &cb[(nl2 + nl3 + nl4)]; // nl33
double *cb34 = &cb[(nl2 + nl3 + nl4 + nl33)]; // nl34
double *cb44 = &cb[(nl2 + nl3 + nl4 + nl33 + nl34)]; // nl44
if ((nl33>0) && (Nj>3)) {
crossdesc_reduction(cb3, cb3, cb33, d3, d3, ind33l, ind33r, nl33);
crossdesc_reduction(cb3, cb3, cb33, d3, d3, ind33l, ind33r, nl33);
}
if ((nl34>0) && (Nj>4)) {
crossdesc_reduction(cb3, cb4, cb34, d3, d4, ind34l, ind34r, nl34);
}
crossdesc_reduction(cb3, cb4, cb34, d3, d4, ind34l, ind34r, nl34);
}
if ((nl44>0) && (Nj>5)) {
crossdesc_reduction(cb4, cb4, cb44, d4, d4, ind44l, ind44r, nl44);
}
if ((nl2 > 0) && (Nj>0)) twobody_forces(fij, cb2, rbfx, rbfy, rbfz, tj, Nj);
crossdesc_reduction(cb4, cb4, cb44, d4, d4, ind44l, ind44r, nl44);
}
if ((nl2 > 0) && (Nj>0)) twobody_forces(fij, cb2, rbfx, rbfy, rbfz, tj, Nj);
// print_matrix("cb2", 1, nl2, cb2, 1);
// print_matrix("rbfx", Nj, nrbf2, rbfx, Nj);
// print_matrix("rbfy", Nj, nrbf2, rbfy, Nj);
// print_matrix("rbfz", Nj, nrbf2, rbfz, Nj);
// print_matrix("fij", 3, Nj, fij, 3);
// error->all(FLERR,"stop");
// Initialize forcecoeff to zero
// Initialize forcecoeff to zero
double *forcecoeff = &cb[(nl2 + nl3 + nl4)]; // nl33
std::fill(forcecoeff, forcecoeff + nelements * K3 * nrbf3, 0.0);
if ((nl3 > 0) && (Nj>1)) threebody_forcecoeff(forcecoeff, cb3, sumU);
if ((nl4 > 0) && (Nj>2)) fourbody_forcecoeff(forcecoeff, cb4, sumU);
if ((nl3 > 0) && (Nj>1)) allbody_forces(fij, forcecoeff, Ux, Uy, Uz, tj, Nj);
if ((nl3 > 0) && (Nj>1)) allbody_forces(fij, forcecoeff, Ux, Uy, Uz, tj, Nj);
return e;
}
@ -1222,51 +1222,51 @@ double EAPOD::peratomenergyforce(double *fij, double *rij, double *temp,
int *ti, int *tj, int Nj)
{
if (Nj==0) {
return coeff[nCoeffPerElement*(ti[0]-1)];
return coeff[nCoeffPerElement*(ti[0]-1)];
}
int N = 3*Nj;
for (int n=0; n<N; n++) fij[n] = 0.0;
double *coeff1 = &coeff[nCoeffPerElement*(ti[0]-1)];
double e = coeff1[0];
// calculate base descriptors and their derivatives with respect to atom coordinates
peratombase_descriptors(bd, bdd, rij, temp, ti, tj, Nj);
peratombase_descriptors(bd, bdd, rij, temp, ti, tj, Nj);
if (nClusters > 1) { // multi-environment descriptors
// calculate multi-environment descriptors and their derivatives with respect to atom coordinates
peratomenvironment_descriptors(pd, pdd, bd, bdd, temp, ti[0] - 1, Nj);
for (int j = 0; j<nClusters; j++)
for (int m=0; m<Mdesc; m++)
for (int m=0; m<Mdesc; m++)
e += coeff1[1 + m + j*Mdesc]*bd[m]*pd[j];
double *cb = &temp[0];
double *cp = &temp[Mdesc];
for (int m = 0; m<Mdesc; m++) cb[m] = 0.0;
for (int m = 0; m<Mdesc; m++) cb[m] = 0.0;
for (int j = 0; j<nClusters; j++) cp[j] = 0.0;
for (int j = 0; j<nClusters; j++)
for (int m = 0; m<Mdesc; m++) {
for (int j = 0; j<nClusters; j++)
for (int m = 0; m<Mdesc; m++) {
cb[m] += coeff1[1 + m + j*Mdesc]*pd[j];
cp[j] += coeff1[1 + m + j*Mdesc]*bd[m];
}
char chn = 'N';
char chn = 'N';
double alpha = 1.0, beta = 0.0;
int inc1 = 1;
DGEMV(&chn, &N, &Mdesc, &alpha, bdd, &N, cb, &inc1, &beta, fij, &inc1);
DGEMV(&chn, &N, &nClusters, &alpha, pdd, &N, cp, &inc1, &alpha, fij, &inc1);
DGEMV(&chn, &N, &Mdesc, &alpha, bdd, &N, cb, &inc1, &beta, fij, &inc1);
DGEMV(&chn, &N, &nClusters, &alpha, pdd, &N, cp, &inc1, &alpha, fij, &inc1);
}
else { // single-environment descriptors
for (int m=0; m<Mdesc; m++)
for (int m=0; m<Mdesc; m++)
e += coeff1[1+m]*bd[m];
char chn = 'N';
char chn = 'N';
double alpha = 1.0, beta = 0.0;
int inc1 = 1;
DGEMV(&chn, &N, &Mdesc, &alpha, bdd, &N, &coeff1[1], &inc1, &beta, fij, &inc1);
DGEMV(&chn, &N, &Mdesc, &alpha, bdd, &N, &coeff1[1], &inc1, &beta, fij, &inc1);
}
return e;
}
@ -1333,7 +1333,7 @@ void EAPOD::base_descriptors(double *basedesc, double *x,
int *tj = &tmpint[3*Nj]; // Nj
myneighbors(rij, x, ai, aj, ti, tj, jlist, pairnumsum, atomtype, alist, i);
// many-body descriptors
peratombase_descriptors(bd, bdd, rij, &tmpmem[3*Nj], ti, tj, Nj);
@ -1342,7 +1342,7 @@ void EAPOD::base_descriptors(double *basedesc, double *x,
}
}
}
}
}
void EAPOD::descriptors(double *gd, double *gdd, double *basedesc, double *x,
@ -1376,7 +1376,7 @@ void EAPOD::descriptors(double *gd, double *gdd, double *basedesc, double *x,
int *tj = &tmpint[3*Nj]; // Nj
myneighbors(rij, x, ai, aj, ti, tj, jlist, pairnumsum, atomtype, alist, i);
// many-body descriptors
peratombase_descriptors(bd, bdd, rij, &tmpmem[3*Nj], ti, tj, Nj);
@ -1398,7 +1398,7 @@ void EAPOD::descriptors(double *gd, double *gdd, double *basedesc, double *x,
}
}
}
}
}
void EAPOD::descriptors(double *gd, double *gdd, double *basedesc, double *probdesc, double *x,
@ -1433,7 +1433,7 @@ void EAPOD::descriptors(double *gd, double *gdd, double *basedesc, double *probd
int *tj = &tmpint[3*Nj]; // Nj
myneighbors(rij, x, ai, aj, ti, tj, jlist, pairnumsum, atomtype, alist, i);
// many-body descriptors
peratombase_descriptors(bd, bdd, rij, &tmpmem[3*Nj], ti, tj, Nj);
@ -1462,7 +1462,7 @@ void EAPOD::descriptors(double *gd, double *gdd, double *basedesc, double *probd
}
}
}
}
}
void EAPOD::fourbodydesc23(double *d23, double *d2, double *d3)
@ -1494,15 +1494,15 @@ void EAPOD::crossdescderiv(double *dd12, double *d1, double *d2, double *dd1, do
dd12[n + N*i] = d1[ind1[i]]*dd2[n + N*ind2[i]] + dd1[n + N*ind1[i]]*d2[ind2[i]];
}
void EAPOD::crossdesc_reduction(double *cb1, double *cb2, double *c12, double *d1,
void EAPOD::crossdesc_reduction(double *cb1, double *cb2, double *c12, double *d1,
double *d2, int *ind1, int *ind2, int n12)
{
{
for (int m = 0; m < n12; m++) {
int k1 = ind1[m]; // dd1
int k2 = ind2[m]; // dd2
double c = c12[m];
cb1[k1] += c * d2[k2];
cb2[k2] += c * d1[k1];
double c = c12[m];
cb1[k1] += c * d2[k2];
cb2[k2] += c * d1[k1];
}
}
@ -1550,14 +1550,14 @@ void EAPOD::fourbodydesc(double *d4, double *sumU)
for (int i3 = i2; i3 < nelements; i3++) {
double c3 = sumU[idxU + i3 + nelements * j3];
int kk = p + nabf4 * m + nabf4 * nrbf4 * k;
d4[kk] += t12 * c3;
d4[kk] += t12 * c3;
k += 1;
}
}
}
}
}
}
}
}
void EAPOD::fourbodydescderiv(double *d4, double *dd4, double *sumU, double *Ux, double *Uy,
@ -1587,7 +1587,7 @@ void EAPOD::fourbodydescderiv(double *d4, double *dd4, double *sumU, double *Ux,
double c2 = c*sumU[j2 + K4*m];
double t12 = c1*sumU[j2 + K4*m];
double c3 = sumU[j3 + K4*m];
double t13 = c1*c3;
double t13 = c1*c3;
double t23 = c2*c3;
int kk = p + nabf4*m;
int ii = 3*N*(p + nabf4*m);
@ -1956,7 +1956,7 @@ void EAPOD::radialbasis(double *rbf, double *rbfx, double *rbfy, double *rbfz, d
for (int i=0; i<inversedegree; i++) {
int p = besseldegree*nbesselpars + i;
int nij = n + N*p;
double a = pow(dij, (double) (i+1.0));
double a = powint(dij, i+1);
rbf[nij] = fcut/a;
@ -2219,7 +2219,7 @@ void EAPOD::snapshots(double *rbf, double *xij, int N)
double y2 = y*y;
double y3 = 1.0 - y2*y;
double y4 = y3*y3 + 1e-6;
double y5 = pow(y4, 0.5);
double y5 = sqrt(y4);
double y6 = exp(-1.0/y5);
// Compute the cutoff function
@ -2493,7 +2493,7 @@ void EAPOD::allocate_temp_memory(int Nj)
memory->create(bd, Mdesc, "bdd");
memory->create(bdd, 3*Nj*Mdesc, "bdd");
memory->create(pd, nClusters, "bdd");
memory->create(pdd, 3*Nj*nClusters, "bdd");
memory->create(pdd, 3*Nj*nClusters, "bdd");
}
void EAPOD::free_temp_memory()
@ -2684,7 +2684,7 @@ void EAPOD::MatMul(double *c, double *a, double *b, int r1, int c1, int c2)
c[i + r1*j] += a[i + r1*k] * b[k + c1*j];
}
void EAPOD::peratomenvironment_descriptors(double *P, double *dP_dR, double *B, double *dB_dR, double *tmp, int elem, int nNeighbors)
void EAPOD::peratomenvironment_descriptors(double *P, double *dP_dR, double *B, double *dB_dR, double *tmp, int elem, int nNeighbors)
{
double *ProjMat = &Proj[nComponents*Mdesc*elem];
double *centroids = &Centroids[nComponents*nClusters*elem];
@ -2693,7 +2693,7 @@ void EAPOD::peratomenvironment_descriptors(double *P, double *dP_dR, double *B,
double *dD_dpca = &tmp[nComponents + nClusters];
double *dD_dB = &tmp[nComponents + nClusters + nClusters*nComponents];
double *dP_dD = &tmp[nComponents + nClusters + nClusters*nComponents + nClusters*Mdesc];
double *dP_dB = &tmp[nComponents + nClusters + nClusters*nComponents + nClusters*Mdesc + nClusters*nClusters];
double *dP_dB = &tmp[nComponents + nClusters + nClusters*nComponents + nClusters*Mdesc + nClusters*nClusters];
// calculate principal components
for (int k = 0; k < nComponents; k++) {
@ -2731,7 +2731,7 @@ void EAPOD::peratomenvironment_descriptors(double *P, double *dP_dR, double *B,
char cht = 'T';
double alpha = 1.0, beta = 0.0;
DGEMM(&chn, &chn, &nClusters, &Mdesc, &nComponents, &alpha, dD_dpca, &nClusters, ProjMat, &nComponents, &beta, dD_dB, &nClusters);
// calculate dP_dD
double S1 = 1 / sumD;
double S2 = sumD * sumD;
@ -2744,12 +2744,12 @@ void EAPOD::peratomenvironment_descriptors(double *P, double *dP_dR, double *B,
dP_dD[j + j * nClusters] += S1;
}
// calculate dP_dB = dP_dD * dD_dB, which are derivatives of probabilities with respect to local descriptors
// calculate dP_dB = dP_dD * dD_dB, which are derivatives of probabilities with respect to local descriptors
DGEMM(&chn, &chn, &nClusters, &Mdesc, &nClusters, &alpha, dP_dD, &nClusters, dD_dB, &nClusters, &beta, dP_dB, &nClusters);
// calculate dP_dR = dB_dR * dP_dB , which are derivatives of probabilities with respect to atomic coordinates
// calculate dP_dR = dB_dR * dP_dB , which are derivatives of probabilities with respect to atomic coordinates
int N = 3*nNeighbors;
DGEMM(&chn, &cht, &N, &nClusters, &Mdesc, &alpha, dB_dR, &N, dP_dB, &nClusters, &beta, dP_dR, &N);
DGEMM(&chn, &cht, &N, &nClusters, &Mdesc, &alpha, dB_dR, &N, dP_dB, &nClusters, &beta, dP_dR, &N);
}
void EAPOD::init3bodyarray(int *np, int *pq, int *pc, int Pa)

26
src/ML-POD/eapod.h
View File

@ -86,7 +86,7 @@ public:
double rin;
double rcut;
int true4BodyDesc;
int nelements; // number of elements
int pbc[3];
@ -110,8 +110,8 @@ public:
int nClusters; // number of environment clusters
int nComponents; // number of principal components
//int nNeighbors; // numbe of neighbors
int Mdesc; // number of base descriptors
int Mdesc; // number of base descriptors
double *Proj; // PCA Projection matrix
double *Centroids; // centroids of the clusters
double *bd; // base descriptors
@ -125,7 +125,7 @@ public:
int Njmax;
int nCoeffPerElement; // number of coefficients per element = (nl1 + Mdesc*nClusters)
int nCoeffAll; // number of coefficients for all elements = (nl1 + Mdesc*nClusters)*nelements
int ncoeff; // number of coefficients in the input file
int ncoeff; // number of coefficients in the input file
int ns; // number of snapshots for radial basis functions
int nd1, nd2, nd3, nd4, nd5, nd6, nd7, nd; // number of global descriptors
int nl1, nl2, nl3, nl4, nl5, nl6, nl7, nl; // number of local descriptors
@ -167,19 +167,19 @@ public:
int read_projection_matrix(std::string proj_file);
int read_centroids(std::string centroids_file);
int estimate_temp_memory(int Nj);
void free_temp_memory();
int estimate_temp_memory(int Nj);
void free_temp_memory();
void allocate_temp_memory(int Nj);
//void mknewcoeff();
void mknewcoeff(double *c, int nc);
void twobodydesc(double *d2, double *rbf, int *tj, int N);
void twobodydesc(double *d2, double *rbf, int *tj, int N);
void twobodydescderiv(double *d2, double *dd2, double *rbf, double *rbfx,
double *rbfy, double *rbfz, int *tj, int N);
void twobody_forces(double *fij, double *cb2, double *rbfx, double *rbfy, double *rbfz, int *tj, int Nj);
void threebodydesc(double *d3, double *sumU);
void threebodydescderiv(double *dd3, double *sumU, double *Ux, double *Uy, double *Uz,
int *atomtype, int N);
@ -189,7 +189,7 @@ public:
void fourbodydescderiv(double *d4, double *dd4, double *sumU, double *Ux, double *Uy, double *Uz,
int *atomtype, int N);
void fourbody_forcecoeff(double *fb4, double *cb4, double *sumU);
void allbody_forces(double *fij, double *forcecoeff, double *rbf, double *rbfx, double *rbfy, double *rbfz,
double *abf, double *abfx, double *abfy, double *abfz, int *tj, int Nj);
void allbody_forces(double *fij, double *forcecoeff, double *Ux, double *Uy, double *Uz, int *tj, int Nj);
@ -204,7 +204,7 @@ public:
int *ti, int *tj, int Nj);
double peratombase_coefficients(double *cb, double *bd, int *ti);
double peratom_environment_descriptors(double *cb, double *bd, double *tm, int *ti);
void peratomenvironment_descriptors(double *P, double *dP_dR, double *B, double *dB_dR, double *tmp, int elem, int nNeighbors);
void base_descriptors(double *basedesc, double *x, int *atomtype, int *alist,
@ -212,7 +212,7 @@ public:
void descriptors(double *basedesc, double *probdesc, double *x, int *atomtype, int *alist,
int *jlist, int *pairnumsum, int natom);
double peratomenergyforce(double *fij, double *rij, double *temp, int *ti, int *tj, int Nj);
double peratomenergyforce2(double *fij, double *rij, double *temp, int *ti, int *tj, int Nj);
@ -226,8 +226,8 @@ public:
void crossdesc(double *d12, double *d1, double *d2, int *ind1, int *ind2, int n12);
void crossdescderiv(double *dd12, double *d1, double *d2, double *dd1, double *dd2,
int *ind1, int *ind2, int n12, int N);
void crossdesc_reduction(double *cb1, double *cb2, double *c12, double *d1,
int *ind1, int *ind2, int n12, int N);
void crossdesc_reduction(double *cb1, double *cb2, double *c12, double *d1,
double *d2, int *ind1, int *ind2, int n12);
};

126
src/ML-POD/fitpod_command.cpp
View File

@ -64,9 +64,9 @@ void FitPOD::command(int narg, char **arg)
cent_file = std::string(arg[4]); // centroid input file
else
cent_file = "";
fastpodptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, cent_file);
fastpodptr = new EAPOD(lmp, pod_file, coeff_file, proj_file, cent_file);
desc.nCoeffAll = fastpodptr->nCoeffAll;
desc.nClusters = fastpodptr->nClusters;
read_data_files(data_file, fastpodptr->species);
@ -76,11 +76,11 @@ void FitPOD::command(int narg, char **arg)
if (desc.nClusters > 1) estimate_memory_neighborstruct(envdata, fastpodptr->pbc, fastpodptr->rcut, fastpodptr->nelements);
allocate_memory_neighborstruct();
estimate_memory_fastpod(traindata);
estimate_memory_fastpod(testdata);
estimate_memory_fastpod(testdata);
allocate_memory_descriptorstruct(fastpodptr->nCoeffAll);
if (coeff_file != "") podArrayCopy(desc.c, fastpodptr->coeff, fastpodptr->nCoeffAll);
if (compute_descriptors==0) {
if (((int) envdata.data_path.size() > 1) && (desc.nClusters > 1)) {
@ -631,7 +631,7 @@ void FitPOD::get_data(datastruct &data, const std::vector<std::string> &species)
// Group weights have same size as energy.
memory->create(data.we, n, "fitpod:we");
memory->create(data.wf, n, "fitpod:wf");
n = data.num_atom_sum;
memory->create(data.position, 3*n, "fitpod:position");
memory->create(data.force, 3*n, "fitpod:force");
@ -811,7 +811,7 @@ void FitPOD::select_data(datastruct &newdata, const datastruct &data)
// Group weights have same size as energy.
memory->create(newdata.we, n, "fitpod:we");
memory->create(newdata.wf, n, "fitpod:wf");
n = newdata.num_atom_sum;
memory->create(newdata.position, 3*n, "fitpod:newdata_position");
memory->create(newdata.force, 3*n, "fitpod:newdata_force");
@ -942,7 +942,7 @@ void FitPOD::read_data_files(const std::string& data_file, const std::vector<std
compute_descriptors = 1;
if (comm->me == 0)
utils::logmesg(lmp, "**************** Begin of Environment Configuration Set ****************\n");
get_data(envdata, species);
get_data(envdata, species);
if (comm->me == 0)
utils::logmesg(lmp, "**************** End of Environment Configuration Set ****************\n");
compute_descriptors = tmp;
@ -993,7 +993,7 @@ void FitPOD::read_data_files(const std::string& data_file, const std::vector<std
}
else {
testdata.data_path = traindata.data_path;
}
}
}
int FitPOD::latticecoords(double *y, int *alist, double *x, double *a1, double *a2, double *a3, double rcut, int *pbc, int nx)
@ -1121,7 +1121,7 @@ void FitPOD::allocate_memory_neighborstruct()
}
void FitPOD::allocate_memory_descriptorstruct(int nCoeffAll)
{
{
memory->create(desc.bd, nb.natom_max*fastpodptr->Mdesc, "fitpod:desc_ld");
memory->create(desc.pd, nb.natom_max*fastpodptr->nClusters, "fitpod:desc_ld");
memory->create(desc.gd, nCoeffAll, "fitpod:desc_gd");
@ -1280,21 +1280,21 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
{
if (comm->me == 0)
utils::logmesg(lmp, "**************** Begin Calculating Environment Descriptor Matrix ****************\n");
//printf("number of configurations = %d\n", (int) data.num_atom.size());
//printf("number of configurations = %d\n", (int) data.num_atom.size());
int nComponents = fastpodptr->nComponents;
int Mdesc = fastpodptr->Mdesc;
int nClusters = fastpodptr->nClusters;
int nelements = fastpodptr->nelements;
int nelements = fastpodptr->nelements;
memory->create(fastpodptr->Centroids, nClusters*nComponents*nelements, "fitpod:centroids");
memory->create(fastpodptr->Proj, Mdesc*nComponents*nelements, "fitpod:P");
int nAtoms = 0;
int nTotalAtoms = 0;
int nTotalAtoms = 0;
for (int ci=0; ci < (int) data.num_atom.size(); ci++) {
if ((ci % comm->nprocs) == comm->me) nAtoms += data.num_atom[ci];
nTotalAtoms += data.num_atom[ci];
if ((ci % comm->nprocs) == comm->me) nAtoms += data.num_atom[ci];
nTotalAtoms += data.num_atom[ci];
}
double *basedescmatrix = (double *) malloc(nAtoms*Mdesc*sizeof(double));
@ -1314,19 +1314,19 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
char chu = 'U';
double alpha = 1.0, beta = 0.0;
for (int elem=0; elem < nelements; elem++) {
for (int elem=0; elem < nelements; elem++) {
nElemAtoms[elem] = 0; // number of atoms for this element
}
for (int ci=0; ci < (int) data.num_atom.size(); ci++) {
if ((ci % comm->nprocs) == comm->me) {
if ((ci % comm->nprocs) == comm->me) {
int natom = data.num_atom[ci];
int natom_cumsum = data.num_atom_cumsum[ci];
int *atomtype = &data.atomtype[natom_cumsum];
for (int n=0; n<natom; n++)
nElemAtoms[atomtype[n]-1] += 1;
for (int n=0; n<natom; n++)
nElemAtoms[atomtype[n]-1] += 1;
}
}
nElemAtomsCumSum[0] = 0;
for (int elem=0; elem < nelements; elem++) {
nElemAtomsCumSum[elem+1] = nElemAtomsCumSum[elem] + nElemAtoms[elem];
@ -1338,7 +1338,7 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
// printf("%d %d %d\n",comm->me, nElemAtoms[0], nElemAtoms[1]);
// printf("%d %d %d\n",comm->me, nElemAtomsCount[0], nElemAtomsCount[1]);
// printf("%d %d %d %d\n",comm->me, nElemAtomsCumSum[0], nElemAtomsCumSum[1], nElemAtomsCumSum[2]);
MPI_Barrier(MPI_COMM_WORLD);
// loop over each configuration in the data set
@ -1352,45 +1352,45 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
base_descriptors_fastpod(data, ci);
// basedescmatrix is a Mdesc x nAtoms matrix
int natom = data.num_atom[ci];
int natom = data.num_atom[ci];
int natom_cumsum = data.num_atom_cumsum[ci];
int *atomtype = &data.atomtype[natom_cumsum];
for (int n=0; n<natom; n++) {
int elem = atomtype[n]-1; // offset by 1 to match the element index in the C++ code
int elem = atomtype[n]-1; // offset by 1 to match the element index in the C++ code
nElemAtomsCount[elem] += 1;
int k = nElemAtomsCumSum[elem] + nElemAtomsCount[elem] - 1;
for (int m=0; m<Mdesc; m++)
basedescmatrix[m + Mdesc*k] = desc.bd[n + natom*(m)];
}
for (int m=0; m<Mdesc; m++)
basedescmatrix[m + Mdesc*k] = desc.bd[n + natom*(m)];
}
}
}
MPI_Barrier(MPI_COMM_WORLD);
// printf("%d %d %d\n",comm->me, nElemAtomsCount[0], nElemAtomsCount[1]);
// MPI_Barrier(MPI_COMM_WORLD);
//
// for (int elem=0; elem < nelements; elem++) { // loop over each element
//
// for (int elem=0; elem < nelements; elem++) { // loop over each element
// nAtoms = nElemAtoms[elem];
// nTotalAtoms = nAtoms;
//
//
// MPI_Barrier(MPI_COMM_WORLD);
// MPI_Allreduce(MPI_IN_PLACE, &nTotalAtoms, 1, MPI_INT, MPI_SUM, world);
//
//
// double *descmatrix = &basedescmatrix[Mdesc*nElemAtomsCumSum[elem]];
// DGEMM(&chn, &cht, &Mdesc, &Mdesc, &nAtoms, &alpha, descmatrix, &Mdesc, descmatrix, &Mdesc, &beta, A, &Mdesc);
//
// MPI_Barrier(MPI_COMM_WORLD);
// MPI_Allreduce(MPI_IN_PLACE, A, Mdesc*Mdesc, MPI_DOUBLE, MPI_SUM, world);
//
//
// MPI_Barrier(MPI_COMM_WORLD);
// MPI_Allreduce(MPI_IN_PLACE, A, Mdesc*Mdesc, MPI_DOUBLE, MPI_SUM, world);
//
// if (comm->me == 0) print_matrix("A", Mdesc, Mdesc, A, Mdesc);
// }
//
//
// MPI_Barrier(MPI_COMM_WORLD);
// error->all(FLERR, "stop for debugging");
int save = 0;
for (int elem=0; elem < nelements; elem++) { // loop over each element
for (int elem=0; elem < nelements; elem++) { // loop over each element
nAtoms = nElemAtoms[elem];
nTotalAtoms = nAtoms;
MPI_Allreduce(MPI_IN_PLACE, &nTotalAtoms, 1, MPI_INT, MPI_SUM, world);
@ -1400,11 +1400,11 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
double *centroids = &fastpodptr->Centroids[nComponents*nClusters*elem];
// Calculate covariance matrix A = basedescmatrix*basedescmatrix'. A is a Mdesc x Mdesc matrix
DGEMM(&chn, &cht, &Mdesc, &Mdesc, &nAtoms, &alpha, descmatrix, &Mdesc, descmatrix, &Mdesc, &beta, A, &Mdesc);
MPI_Allreduce(MPI_IN_PLACE, A, Mdesc*Mdesc, MPI_DOUBLE, MPI_SUM, world);
DGEMM(&chn, &cht, &Mdesc, &Mdesc, &nAtoms, &alpha, descmatrix, &Mdesc, descmatrix, &Mdesc, &beta, A, &Mdesc);
MPI_Allreduce(MPI_IN_PLACE, A, Mdesc*Mdesc, MPI_DOUBLE, MPI_SUM, world);
//if (comm->me == 0) print_matrix("A", Mdesc, Mdesc, A, Mdesc);
if ((comm->me == 0) && (save==1))
savematrix2binfile(data.filenametag + "_covariance_matrix_elem" + std::to_string(elem+1) + ".bin", A, Mdesc, Mdesc);
@ -1427,15 +1427,15 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
// Calculate principal compoment analysis matrix pca = P*descmatrix. pca is a nComponents x nAtoms matrix
DGEMM(&chn, &chn, &nComponents, &nAtoms, &Mdesc, &alpha, Proj, &nComponents, descmatrix, &Mdesc, &beta, pca, &nComponents);
// initialize centroids
for (int i = 0; i < nClusters * nComponents; i++) centroids[i] = 0.0;
for (int i=0; i < nAtoms; i++) {
// initialize centroids
for (int i = 0; i < nClusters * nComponents; i++) centroids[i] = 0.0;
for (int i=0; i < nAtoms; i++) {
int m = (i*nClusters)/nAtoms;
for (int j=0; j < nComponents; j++)
centroids[j + nComponents*m] += pca[j + nComponents*i];
}
MPI_Allreduce(MPI_IN_PLACE, centroids, nClusters * nComponents, MPI_DOUBLE, MPI_SUM, world);
MPI_Allreduce(MPI_IN_PLACE, centroids, nClusters * nComponents, MPI_DOUBLE, MPI_SUM, world);
double fac = ((double) nClusters)/((double) nTotalAtoms);
for (int i = 0; i < nClusters * nComponents; i++) centroids[i] = centroids[i]*fac;
//for (int i = 0; i < desc.nClusters * nComponents; i++) printf("centroids[%d] = %f\n", i, centroids[i]);
@ -1449,21 +1449,21 @@ void FitPOD::environment_cluster_calculation(const datastruct &data)
savematrix2binfile(data.filenametag + "_eigenvector_matrix_elem" + std::to_string(elem+1) + ".bin", A, Mdesc, Mdesc);
savematrix2binfile(data.filenametag + "_eigenvalues_elem" + std::to_string(elem+1) + ".bin", b, Mdesc, 1);
}
savematrix2binfile(data.filenametag + "_desc_matrix_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", descmatrix, Mdesc, nAtoms);
savematrix2binfile(data.filenametag + "_desc_matrix_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", descmatrix, Mdesc, nAtoms);
savematrix2binfile(data.filenametag + "_pca_matrix_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", pca, nComponents, nAtoms);
saveintmatrix2binfile(data.filenametag + "_cluster_assignments_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", assignments, nAtoms, 1);
}
//savematrix2binfile(data.filenametag + "_pca_matrix_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", pca, nComponents, nAtoms);
//saveintmatrix2binfile(data.filenametag + "_cluster_assignments_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", assignments, nAtoms, 1);
//saveintmatrix2binfile(data.filenametag + "_cluster_assignments_elem" + std::to_string(elem+1) + "_proc" + std::to_string(comm->me+1) + ".bin", assignments, nAtoms, 1);
MPI_Barrier(MPI_COMM_WORLD);
}
if (comm->me == 0) {
savedata2textfile(data.filenametag + "_projection_matrix" + ".pod", "projection_matrix: {}\n ", fastpodptr->Proj, nComponents*Mdesc*nelements, 1, 1);
savedata2textfile(data.filenametag + "_centroids" + ".pod", "centroids: {} \n", fastpodptr->Centroids, nComponents*nClusters*nelements, 1, 1);
savedata2textfile(data.filenametag + "_centroids" + ".pod", "centroids: {} \n", fastpodptr->Centroids, nComponents*nClusters*nelements, 1, 1);
}
free(basedescmatrix);
free(pca);
free(A);
@ -1827,9 +1827,9 @@ void FitPOD::error_analysis(const datastruct &data, double *coeff)
int natom = data.num_atom[ci];
int nforce = dim*natom;
emae += outarray[4 + m*ci]; // sum_c |ePOD_c - eDFT_c|/natom_c
essr += outarray[4 + m*ci]*outarray[4 + m*ci]; // sum_c |ePOD_c - eDFT_c|^2/(natom_c)^2
fmae += outarray[7 + m*ci]*nforce; // sum_c |fPOD_c - fDFT_c|
emae += outarray[4 + m*ci]; // sum_c |ePOD_c - eDFT_c|/natom_c
essr += outarray[4 + m*ci]*outarray[4 + m*ci]; // sum_c |ePOD_c - eDFT_c|^2/(natom_c)^2
fmae += outarray[7 + m*ci]*nforce; // sum_c |fPOD_c - fDFT_c|
fssr += ssrarray[ci];
nforceall += nforce;
ci += 1;
@ -1844,15 +1844,15 @@ void FitPOD::error_analysis(const datastruct &data, double *coeff)
errors[3 + 4*q] = sqrt(fssr/nforceall);
nf += nforceall;
errors[0] += emae; // sum_c |ePOD_c - eDFT_c|/natom_c
errors[1] += essr; // sum_c |ePOD_c - eDFT_c|^2/(natom_c)^2
errors[0] += emae; // sum_c |ePOD_c - eDFT_c|/natom_c
errors[1] += essr; // sum_c |ePOD_c - eDFT_c|^2/(natom_c)^2
errors[2] += fmae;
errors[3] += fssr;
}
if (nc == 0) nc = 1;
if (nf == 0) nf = 1;
errors[0] = errors[0]/nc; // (1/Nc) * sum_c |ePOD_c - eDFT_c|/natom_c
errors[0] = errors[0]/nc; // (1/Nc) * sum_c |ePOD_c - eDFT_c|/natom_c
errors[1] = sqrt(errors[1]/nc); // sqrt { (1/Nc) * sum_c |ePOD_c - eDFT_c|^2/(natom_c)^2 }
errors[2] = errors[2]/nf;
errors[3] = sqrt(errors[3]/nf);
@ -2210,13 +2210,13 @@ void FitPOD::saveintmatrix2binfile(std::string filename, int *A, int nrows, int
void FitPOD::savedata2textfile(std::string filename, std::string text, double *A, int n, int m, int dim)
{
if (comm->me == 0) {
if (comm->me == 0) {
int precision = 15;
FILE *fp = fopen(filename.c_str(), "w");
FILE *fp = fopen(filename.c_str(), "w");
if (dim==1) {
fmt::print(fp, text, n);
for (int i = 0; i < n; i++)
fmt::print(fp, "{:<10.{}f} \n", A[i], precision);
fmt::print(fp, "{:<10.{}f} \n", A[i], precision);
}
else if (dim==2) {
fmt::print(fp, text, n);

14
src/ML-POD/fitpod_command.h
View File

@ -124,14 +124,14 @@ private:
double *gdd=nullptr; // derivatives of global descriptors and peratom descriptors
double *A=nullptr; // least-square matrix for all descriptors
double *b=nullptr; // least-square vector for all descriptors
double *c=nullptr; // coefficents of descriptors
int szd = 0;
double *c=nullptr; // coefficents of descriptors
int szd = 0;
int nCoeffAll = 0; // number of global descriptors
int nClusters = 0; // number of environment clusters
int nClusters = 0; // number of environment clusters
};
int save_descriptors = 0;
int compute_descriptors = 0;
int compute_descriptors = 0;
datastruct traindata;
datastruct testdata;
datastruct envdata;
@ -162,7 +162,7 @@ private:
void matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom);
void matrix33_inverse(double *invA, double *A1, double *A2, double *A3);
double squareDistance(const double *a, const double *b, int DIMENSIONS);
double squareDistance(const double *a, const double *b, int DIMENSIONS);
void assignPointsToClusters(double *points, double *centroids, int *assignments, int *clusterSizes, int NUM_POINTS, int NUM_CLUSTERS, int DIMENSION);
void updateCentroids(double *points, double *centroids, int *assignments, int *clusterSizes, int NUM_POINTS, int NUM_CLUSTERS, int DIMENSIONS);
void KmeansClustering(double *points, double *centroids, int *assignments, int *clusterSizes, int NUM_POINTS, int NUM_CLUSTERS, int DIMENSIONS, int MAX_ITER);
@ -196,7 +196,7 @@ private:
void estimate_memory_neighborstruct(const datastruct &data, int *pbc, double rcut, int nelements);
void allocate_memory_neighborstruct();
void allocate_memory_descriptorstruct(int nd);
void estimate_memory_fastpod(const datastruct &data);
void estimate_memory_fastpod(const datastruct &data);
void local_descriptors_fastpod(const datastruct &data, int ci);
void base_descriptors_fastpod(const datastruct &data, int ci);
void least_squares_matrix(const datastruct &data, int ci);
@ -204,7 +204,7 @@ private:
void descriptors_calculation(const datastruct &data);
void environment_cluster_calculation(const datastruct &data);
void print_analysis(const datastruct &data, double *outarray, double *errors);
void error_analysis(const datastruct &data, double *coeff);
void error_analysis(const datastruct &data, double *coeff);
double energyforce_calculation_fastpod(double *force, const datastruct &data, int ci);
void energyforce_calculation(const datastruct &data, double *coeff);
};

688
src/ML-POD/pair_pod.cpp
View File

File diff suppressed because it is too large Load Diff

60
src/ML-POD/pair_pod.h
View File

@ -41,13 +41,13 @@ public:
void NeighborCount(double **x, int **firstneigh, int *ilist, int *numneigh, double rcutsq, int i1, int i2);
void NeighborList(double **x, int **firstneigh, int *atomtype, int *map, int *ilist, int *numneigh,
double rcutsq, int i1, int i2);
void tallyenergy(double *ei, int istart, int Ni);
void tallystress(double *fij, double *rij, int *ai, int *aj, int nlocal, int N);
void tallyenergy(double *ei, int istart, int Ni);
void tallystress(double *fij, double *rij, int *ai, int *aj, int nlocal, int N);
void tallyforce(double **force, double *fij, int *ai, int *aj, int N);
void divideInterval(int *intervals, int N, int M);
int calculateNumberOfIntervals(int N, int intervalSize);
int calculateNumberOfIntervals(int N, int intervalSize);
int numberOfNeighbors();
void copy_data_from_pod_class();
void radialbasis(double *rbft, double *rbftx, double *rbfty, double *rbftz, double *rij, int Nij);
void orthogonalradialbasis(int Nij);
@ -65,32 +65,32 @@ public:
void crossdesc(double *d12, double *d1, double *d2, int *ind1, int *ind2, int n12, int Ni);
void crossdescderiv(double *dd12, double *d1, double *d2, double *dd1, double *dd2,
int *ind1, int *ind2, int *idxi, int n12, int Ni, int Nij);
void blockatombase_descriptors(double *bd1, double *bdd1, int Ni, int Nij);
void blockatombase_descriptors(double *bd1, double *bdd1, int Ni, int Nij);
void blockatomenergyforce(double *ei, double *fij, int Ni, int Nij);
void crossdesc_reduction(double *cb1, double *cb2, double *c12, double *d1,
double *d2, int *ind1, int *ind2, int n12, int Ni);
void crossdesc_reduction(double *cb1, double *cb2, double *c12, double *d1,
double *d2, int *ind1, int *ind2, int n12, int Ni);
void blockatom_base_descriptors(double *bd1, int Ni, int Nij);
void blockatom_base_coefficients(double *ei, double *cb, double *B, int Ni);
void blockatom_environment_descriptors(double *ei, double *cb, double *B, int Ni);
void blockatom_energyforce(double *ei, double *fij, int Ni, int Nij);
void blockatom_energies(double *ei, int Ni, int Nij);
void blockatom_forces(double *fij, int Ni, int Nij);
void blockatom_forces(double *fij, int Ni, int Nij);
void twobody_forces(double *fij, double *cb2, int Ni, int Nij);
void threebody_forces(double *fij, double *cb3, int Ni, int Nij);
void fourbody_forces(double *fij, double *cb4, int Ni, int Nij);
void threebody_forcecoeff(double *fb3, double *cb3, int Ni);
void fourbody_forcecoeff(double *fb4, double *cb4, int Ni);
void allbody_forces(double *fij, double *forcecoeff, int Nij);
void savematrix2binfile(std::string filename, double *A, int nrows, int ncols);
void saveintmatrix2binfile(std::string filename, int *A, int nrows, int ncols);
void savedatafordebugging();
void saveintmatrix2binfile(std::string filename, int *A, int nrows, int ncols);
void savedatafordebugging();
protected:
class EAPOD *fastpodptr;
class EAPOD *fastpodptr;
virtual void allocate();
void grow_atoms(int Ni);
void grow_pairs(int Nij);
@ -98,13 +98,13 @@ protected:
int atomBlockSize; // size of each atom block
int nAtomBlocks; // number of atoms blocks
int atomBlocks[101]; // atom blocks
int ni; // total number of atoms i
int nij; // total number of pairs (i,j)
int nimax; // maximum number of atoms i
int nijmax; // maximum number of pairs (i,j)
int nelements; // number of elements
int nelements; // number of elements
int onebody; // one-body descriptors
int besseldegree; // degree of Bessel functions
int inversedegree; // degree of inverse functions
@ -113,57 +113,57 @@ protected:
int ns; // number of snapshots for radial basis functions
int nl1, nl2, nl3, nl4, nl23, nl33, nl34, nl44, nl; // number of local descriptors
int nrbf2, nrbf3, nrbf4, nrbfmax; // number of radial basis functions
int nabf3, nabf4; // number of angular basis functions
int nabf3, nabf4; // number of angular basis functions
int K3, K4, Q4; // number of monomials
// environmental variables
int nClusters; // number of environment clusters
int nComponents; // number of principal components
int Mdesc; // number of base descriptors
int Mdesc; // number of base descriptors
double rin; // inner cut-off radius
double rcut; // outer cut-off radius
double rmax; // rcut - rin
double *rij; // (xj - xi) for all pairs (I, J)
double *fij; // force for all pairs (I, J)
double *ei; // energy for each atom I
int *typeai; // types of atoms I only
int *numij; // number of pairs (I, J) for each atom I
int *numij; // number of pairs (I, J) for each atom I
int *idxi; // storing linear indices of atom I for all pairs (I, J)
int *ai; // IDs of atoms I for all pairs (I, J)
int *aj; // IDs of atoms J for all pairs (I, J)
int *ti; // types of atoms I for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
int *tj; // types of atoms J for all pairs (I, J)
double besselparams[3];
double *Phi ; // eigenvectors matrix ns x ns
double *rbf; // radial basis functions nij x nrbfmax
double *rbfx; // x-derivatives of radial basis functions nij x nrbfmax
double *rbf; // radial basis functions nij x nrbfmax
double *rbfx; // x-derivatives of radial basis functions nij x nrbfmax
double *rbfy; // y-derivatives of radial basis functions nij x nrbfmax
double *rbfz; // z-derivatives of radial basis functions nij x nrbfmax
double *rbfz; // z-derivatives of radial basis functions nij x nrbfmax
double *abf; // angular basis functions nij x K3
double *abfx; // x-derivatives of angular basis functions nij x K3
double *abfy; // y-derivatives of angular basis functions nij x K3
double *abfy; // y-derivatives of angular basis functions nij x K3
double *abfz; // z-derivatives of angular basis functions nij x K3
double *abftm ; // angular basis functions 4 x K3
double *sumU; // sum of radial basis functions ni x K3 x nrbfmax x nelements
double *forcecoeff; // force coefficients ni x K3 x nrbfmax x nelements
double *Proj; // PCA Projection matrix
double *Centroids; // centroids of the clusters
double *Centroids; // centroids of the clusters
double *bd; // base descriptors ni x Mdesc
double *cb; // force coefficients for base descriptors ni x Mdesc
double *pd; // environment probability descriptors ni x nClusters
double *bdd; // base descriptors derivatives 3 x nij x Mdesc
double *bdd; // base descriptors derivatives 3 x nij x Mdesc
double *pdd; // environment probability descriptors derivatives 3 x nij x nClusters
double *coefficients; // coefficients nCoeffPerElement x nelements
int *pq3, *pn3, *pc3; // arrays to compute 3-body angular basis functions
int *pa4, *pb4, *pc4; // arrays to compute 4-body angular basis functions
int *pa4, *pb4, *pc4; // arrays to compute 4-body angular basis functions
int *ind33l, *ind33r; // nl33
int *ind34l, *ind34r; // nl34
int *ind44l, *ind44r; // nl44
int *elemindex;
bool peratom_warn; // print warning about missing per-atom energies or stresses
};