Remove unused functions
This commit is contained in:
exapde
@ -69,7 +69,7 @@ void FitPOD::command(int narg, char **arg)
|
||||
|
||||
// get POD coefficients from an input file
|
||||
|
||||
if (coeff_file != "") podptr->podArrayCopy(desc.c, podptr->pod.coeff, podptr->pod.nd);
|
||||
if (coeff_file != "") podArrayCopy(desc.c, podptr->pod.coeff, podptr->pod.nd);
|
||||
|
||||
// compute POD coefficients using least-squares method
|
||||
|
||||
@ -507,13 +507,13 @@ void FitPOD::get_data(datastruct &data, std::vector<std::string> species)
|
||||
}
|
||||
|
||||
int len = data.num_atom.size();
|
||||
data.num_atom_min = podptr->podArrayMin(&data.num_atom[0], len);
|
||||
data.num_atom_max = podptr->podArrayMax(&data.num_atom[0], len);
|
||||
data.num_atom_min = podArrayMin(&data.num_atom[0], len);
|
||||
data.num_atom_max = podArrayMax(&data.num_atom[0], len);
|
||||
data.num_atom_cumsum.resize(len+1);
|
||||
podptr->podCumsum(&data.num_atom_cumsum[0], &data.num_atom[0], len+1);
|
||||
podCumsum(&data.num_atom_cumsum[0], &data.num_atom[0], len+1);
|
||||
|
||||
data.num_config_cumsum.resize(nfiles+1);
|
||||
podptr->podCumsum(&data.num_config_cumsum[0], &data.num_config[0], nfiles+1);
|
||||
podCumsum(&data.num_config_cumsum[0], &data.num_config[0], nfiles+1);
|
||||
|
||||
// convert all structures to triclinic system
|
||||
|
||||
@ -529,11 +529,11 @@ void FitPOD::get_data(datastruct &data, std::vector<std::string> species)
|
||||
double *a2 = &lattice[3];
|
||||
double *a3 = &lattice[6];
|
||||
|
||||
podptr->matrix33_inverse(Qmat, a1, a2, a3);
|
||||
podptr->triclinic_lattice_conversion(a1, a2, a3, a1, a2, a3);
|
||||
podptr->matrix33_multiplication(Qmat, lattice, Qmat, DIM);
|
||||
podptr->matrix33_multiplication(x, Qmat, x, natom);
|
||||
podptr->matrix33_multiplication(f, Qmat, f, natom);
|
||||
matrix33_inverse(Qmat, a1, a2, a3);
|
||||
triclinic_lattice_conversion(a1, a2, a3, a1, a2, a3);
|
||||
matrix33_multiplication(Qmat, lattice, Qmat, DIM);
|
||||
matrix33_multiplication(x, Qmat, x, natom);
|
||||
matrix33_multiplication(f, Qmat, f, natom);
|
||||
}
|
||||
|
||||
if (comm->me == 0) {
|
||||
@ -640,12 +640,12 @@ void FitPOD::select_data(datastruct &newdata, datastruct data)
|
||||
newdata.num_atom_each_file[file] = num_atom_sum;
|
||||
}
|
||||
int len = newdata.num_atom.size();
|
||||
newdata.num_atom_min = podptr->podArrayMin(&newdata.num_atom[0], len);
|
||||
newdata.num_atom_max = podptr->podArrayMax(&newdata.num_atom[0], len);
|
||||
newdata.num_atom_min = podArrayMin(&newdata.num_atom[0], len);
|
||||
newdata.num_atom_max = podArrayMax(&newdata.num_atom[0], len);
|
||||
newdata.num_atom_cumsum.resize(len+1);
|
||||
podptr->podCumsum(&newdata.num_atom_cumsum[0], &newdata.num_atom[0], len+1);
|
||||
podCumsum(&newdata.num_atom_cumsum[0], &newdata.num_atom[0], len+1);
|
||||
newdata.num_atom_sum = newdata.num_atom_cumsum[len];
|
||||
podptr->podCumsum(&newdata.num_config_cumsum[0], &newdata.num_config[0], nfiles+1);
|
||||
podCumsum(&newdata.num_config_cumsum[0], &newdata.num_config[0], nfiles+1);
|
||||
newdata.num_config_sum = newdata.num_atom.size();
|
||||
|
||||
int n = data.num_config_sum;
|
||||
@ -867,7 +867,7 @@ int FitPOD::podfullneighborlist(double *y, int *alist, int *neighlist, int *numn
|
||||
|
||||
nn = podneighborlist(neighlist, numneigh, y, rcutsq, nx, N, dim);
|
||||
|
||||
podptr->podCumsum(numneighsum, numneigh, nx+1);
|
||||
podCumsum(numneighsum, numneigh, nx+1);
|
||||
|
||||
return nn;
|
||||
}
|
||||
@ -879,9 +879,9 @@ void FitPOD::allocate_memory(datastruct data)
|
||||
memory->create(desc.A, nd*nd, "fitpod:desc_A");
|
||||
memory->create(desc.b, nd, "fitpod:desc_b");
|
||||
memory->create(desc.c, nd, "fitpod:desc_c");
|
||||
podptr->podArraySetValue(desc.A, 0.0, nd*nd);
|
||||
podptr->podArraySetValue(desc.b, 0.0, nd);
|
||||
podptr->podArraySetValue(desc.c, 0.0, nd);
|
||||
podArraySetValue(desc.A, 0.0, nd*nd);
|
||||
podArraySetValue(desc.b, 0.0, nd);
|
||||
podArraySetValue(desc.c, 0.0, nd);
|
||||
|
||||
int dim = 3;
|
||||
int natom_max = data.num_atom_max;
|
||||
@ -1192,7 +1192,7 @@ void FitPOD::least_squares_matrix(datastruct data, int ci)
|
||||
|
||||
// least-square matrix for all descriptors: A = A + (we*we)*(gd^T * gd)
|
||||
|
||||
podptr->podKron(desc.A, desc.gd, desc.gd, we2, nd, nd);
|
||||
podKron(desc.A, desc.gd, desc.gd, we2, nd, nd);
|
||||
|
||||
// least-square matrix for all descriptors derivatives: A = A + (wf*wf) * (gdd^T * gdd)
|
||||
|
||||
@ -1340,7 +1340,7 @@ double FitPOD::energyforce_calculation(double *force, double *coeff, datastruct
|
||||
podptr->podNeighPairs(rij, nb.y, idxi, ai, aj, ti, tj, nb.pairnum_cumsum, atomtype, nb.pairlist, nb.alist, natom);
|
||||
|
||||
double *effectivecoeff = &tmpmem[3*Nij]; // 3*Nij
|
||||
podptr->podArraySetValue(effectivecoeff, 0.0, nd1234);
|
||||
podArraySetValue(effectivecoeff, 0.0, nd1234);
|
||||
|
||||
double energy = podptr->energyforce_calculation(force, coeff, effectivecoeff, desc.gd, rij,
|
||||
&tmpmem[3*Nij+nd1234], nb.pairnum_cumsum, atomtype, idxi, ai, aj, ti, tj, natom, Nij);
|
||||
@ -1496,8 +1496,8 @@ void FitPOD::error_analysis(datastruct data, double *coeff)
|
||||
outarray[2 + m*ci] = energy;
|
||||
outarray[3 + m*ci] = DFTenergy;
|
||||
outarray[4 + m*ci] = fabs(DFTenergy-energy)/natom;
|
||||
outarray[5 + m*ci] = podptr->podArrayNorm(force.data(), nforce);
|
||||
outarray[6 + m*ci] = podptr->podArrayNorm(DFTforce, nforce);
|
||||
outarray[5 + m*ci] = podArrayNorm(force.data(), nforce);
|
||||
outarray[6 + m*ci] = podArrayNorm(DFTforce, nforce);
|
||||
|
||||
double diff, sum = 0.0, ssr = 0.0;
|
||||
for (int j=0; j<dim*natom; j++) {
|
||||
@ -1604,3 +1604,216 @@ void FitPOD::energyforce_calculation(datastruct data, double *coeff)
|
||||
if (comm->me == 0)
|
||||
utils::logmesg(lmp, "**************** End of Energy/Force Calculation ****************\n");
|
||||
}
|
||||
|
||||
void FitPOD::print_matrix(const char *desc, int m, int n, double **a, int /*lda*/ )
|
||||
{
|
||||
int i, j;
|
||||
printf( "\n %s\n", desc );
|
||||
|
||||
for( i = 0; i < m; i++ )
|
||||
{
|
||||
for( j = 0; j < n; j++ ) printf( " %6.12f", a[j][i] );
|
||||
printf( "\n" );
|
||||
}
|
||||
}
|
||||
|
||||
void FitPOD::print_matrix(const char *desc, int m, int n, double *a, int lda )
|
||||
{
|
||||
int i, j;
|
||||
printf( "\n %s\n", desc );
|
||||
|
||||
for( i = 0; i < m; i++ )
|
||||
{
|
||||
for( j = 0; j < n; j++ ) printf( " %6.12f", a[i+j*lda] );
|
||||
printf( "\n" );
|
||||
}
|
||||
}
|
||||
|
||||
void FitPOD::print_matrix(const char *desc, int m, int n, int *a, int lda)
|
||||
{
|
||||
int i, j;
|
||||
printf( "\n %s\n", desc );
|
||||
|
||||
for( i = 0; i < m; i++ )
|
||||
{
|
||||
for( j = 0; j < n; j++ ) printf( " %d", a[i+j*lda] );
|
||||
printf( "\n" );
|
||||
}
|
||||
}
|
||||
|
||||
void FitPOD::podArrayFill(int* output, int start, int length)
|
||||
{
|
||||
for (int j = 0; j < length; ++j)
|
||||
output[j] = start + j;
|
||||
}
|
||||
|
||||
double FitPOD::podArraySum(double *a, int n)
|
||||
{
|
||||
double e = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
e += a[i];
|
||||
return e;
|
||||
}
|
||||
|
||||
double FitPOD::podArrayMin(double *a, int n)
|
||||
{
|
||||
double b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]<b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
double FitPOD::podArrayMax(double *a, int n)
|
||||
{
|
||||
double b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]>b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
int FitPOD::podArrayMin(int *a, int n)
|
||||
{
|
||||
int b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]<b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
int FitPOD::podArrayMax(int *a, int n)
|
||||
{
|
||||
int b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]>b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
void FitPOD::podKron(double *C, double *A, double *B, double alpha, int M1, int M2)
|
||||
{
|
||||
int M = M1*M2;
|
||||
for (int idx=0; idx<M; idx++)
|
||||
{
|
||||
int ib = idx%M2;
|
||||
int ia = (idx-ib)/M2;
|
||||
C[idx] += alpha*A[ia]*B[ib];
|
||||
}
|
||||
}
|
||||
|
||||
void FitPOD::podCumsum(int* output, int* input, int length)
|
||||
{
|
||||
output[0] = 0;
|
||||
for (int j = 1; j < length; ++j)
|
||||
output[j] = input[j - 1] + output[j - 1];
|
||||
}
|
||||
|
||||
double FitPOD::podArrayNorm(double *a, int n)
|
||||
{
|
||||
double e = a[0]*a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
e += a[i]*a[i];
|
||||
return sqrt(e);
|
||||
}
|
||||
|
||||
double FitPOD::podArrayErrorNorm(double *a, double *b, int n)
|
||||
{
|
||||
double e = (a[0]-b[0])*(a[0]-b[0]);
|
||||
for (int i=1; i<n; i++)
|
||||
e += (a[i]-b[i])*(a[i]-b[i]);
|
||||
return sqrt(e);
|
||||
}
|
||||
|
||||
void FitPOD::podArraySetValue(double *y, double a, int n)
|
||||
{
|
||||
for (int i=0; i<n; i++)
|
||||
y[i] = a;
|
||||
}
|
||||
|
||||
void FitPOD::podArrayCopy(double *y, double *x, int n)
|
||||
{
|
||||
for (int i=0; i<n; i++)
|
||||
y[i] = x[i];
|
||||
}
|
||||
|
||||
void FitPOD::rotation_matrix(double *Rmat, double alpha, double beta, double gamma)
|
||||
{
|
||||
double ca = cos(alpha);
|
||||
double cb = cos(beta);
|
||||
double cg = cos(gamma);
|
||||
double sa = sin(alpha);
|
||||
double sb = sin(beta);
|
||||
double sg = sin(gamma);
|
||||
|
||||
Rmat[0] = ca*cg*cb - sa*sg;
|
||||
Rmat[3] = -ca*cb*sg - sa*cg;
|
||||
Rmat[6] = ca*sb;
|
||||
|
||||
Rmat[1] = sa*cg*cb + ca*sg;
|
||||
Rmat[4] = -sa*cb*sg + ca*cg;
|
||||
Rmat[7] = sa*sb;
|
||||
|
||||
Rmat[2] = -sb*cg;
|
||||
Rmat[5] = sb*sg;
|
||||
Rmat[8] = cb;
|
||||
}
|
||||
|
||||
void FitPOD::matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom)
|
||||
{
|
||||
double x1, x2, x3;
|
||||
for (int i=0; i < natom; i++) {
|
||||
x1 = x[0 + 3*i];
|
||||
x2 = x[1 + 3*i];
|
||||
x3 = x[2 + 3*i];
|
||||
xrot[0 + 3*i] = Rmat[0]*x1 + Rmat[3]*x2 + Rmat[6]*x3;
|
||||
xrot[1 + 3*i] = Rmat[1]*x1 + Rmat[4]*x2 + Rmat[7]*x3;
|
||||
xrot[2 + 3*i] = Rmat[2]*x1 + Rmat[5]*x2 + Rmat[8]*x3;
|
||||
}
|
||||
}
|
||||
|
||||
void FitPOD::matrix33_inverse(double *invA, double *A1, double *A2, double *A3)
|
||||
{
|
||||
double a11 = A1[0];
|
||||
double a21 = A1[1];
|
||||
double a31 = A1[2];
|
||||
double a12 = A2[0];
|
||||
double a22 = A2[1];
|
||||
double a32 = A2[2];
|
||||
double a13 = A3[0];
|
||||
double a23 = A3[1];
|
||||
double a33 = A3[2];
|
||||
double detA = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
|
||||
|
||||
invA[0] = (a22*a33 - a23*a32)/detA;
|
||||
invA[1] = (a23*a31 - a21*a33)/detA;
|
||||
invA[2] = (a21*a32 - a22*a31)/detA;
|
||||
invA[3] = (a13*a32 - a12*a33)/detA;
|
||||
invA[4] = (a11*a33 - a13*a31)/detA;
|
||||
invA[5] = (a12*a31 - a11*a32)/detA;
|
||||
invA[6] = (a12*a23 - a13*a22)/detA;
|
||||
invA[7] = (a13*a21 - a11*a23)/detA;
|
||||
invA[8] = (a11*a22 - a12*a21)/detA;
|
||||
}
|
||||
|
||||
void FitPOD::triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B, double *C)
|
||||
{
|
||||
double Anorm = sqrt(A[0]*A[0] + A[1]*A[1] + A[2]*A[2]);
|
||||
double Bnorm = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
|
||||
double Cnorm = sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
|
||||
|
||||
double Ahat[3];
|
||||
Ahat[0] = A[0]/Anorm; Ahat[1] = A[1]/Anorm; Ahat[2] = A[2]/Anorm;
|
||||
|
||||
double ax = Anorm;
|
||||
double bx = B[0]*Ahat[0] + B[1]*Ahat[1] + B[2]*Ahat[2]; //dot(B,Ahat);
|
||||
double by = sqrt(Bnorm*Bnorm - bx*bx); //sqrt(Bnorm^2 - bx^2);// #norm(cross(Ahat,B));
|
||||
double cx = C[0]*Ahat[0] + C[1]*Ahat[1] + C[2]*Ahat[2]; // dot(C,Ahat);
|
||||
double cy = (B[0]*C[0] + B[1]*C[1] + B[2]*C[2] - bx*cx)/by; // (dot(B, C) - bx*cx)/by;
|
||||
double cz = sqrt(Cnorm*Cnorm - cx*cx - cy*cy); // sqrt(Cnorm^2 - cx^2 - cy^2);
|
||||
|
||||
a[0] = ax; a[1] = 0.0; a[2] = 0.0;
|
||||
b[0] = bx; b[1] = by; b[2] = 0.0;
|
||||
c[0] = cx; c[1] = cy; c[2] = cz;
|
||||
}
|
||||
|
||||
|
||||
@ -122,6 +122,32 @@ public:
|
||||
|
||||
FitPOD(LAMMPS *lmp) : Command(lmp) {}
|
||||
|
||||
|
||||
// functions for collecting/collating arrays
|
||||
|
||||
void print_matrix(const char *desc, int m, int n, int *a, int lda);
|
||||
void print_matrix(const char *desc, int m, int n, double *a, int lda);
|
||||
void print_matrix(const char *desc, int m, int n, double **a, int lda);
|
||||
void podCumsum(int *output, int *input, int length);
|
||||
double podArrayNorm(double *a, int n);
|
||||
double podArrayErrorNorm(double *a, double *b, int n);
|
||||
void podArraySetValue(double *y, double a, int n);
|
||||
void podArrayCopy(double *y, double *x, int n);
|
||||
void podArrayFill(int *output, int start, int length);
|
||||
double podArrayMin(double *a, int n);
|
||||
double podArrayMax(double *a, int n);
|
||||
double podArraySum(double *a, int n);
|
||||
int podArrayMin(int *a, int n);
|
||||
int podArrayMax(int *a, int n);
|
||||
void podKron(double *C, double *A, double *B, double alpha, int M1, int M2);
|
||||
void rotation_matrix(double *Rmat, double alpha, double beta, double gamma);
|
||||
void triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B,
|
||||
double *C);
|
||||
void matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom);
|
||||
void matrix33_inverse(double *invA, double *A1, double *A2, double *A3);
|
||||
|
||||
// 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);
|
||||
|
||||
@ -108,43 +108,6 @@ MLPOD::~MLPOD()
|
||||
|
||||
// clang-format off
|
||||
|
||||
|
||||
void MLPOD::print_matrix(const char *desc, int m, int n, double **a, int /*lda*/ )
|
||||
{
|
||||
int i, j;
|
||||
printf( "\n %s\n", desc );
|
||||
|
||||
for( i = 0; i < m; i++ )
|
||||
{
|
||||
for( j = 0; j < n; j++ ) printf( " %6.12f", a[j][i] );
|
||||
printf( "\n" );
|
||||
}
|
||||
}
|
||||
|
||||
void MLPOD::print_matrix(const char *desc, int m, int n, double *a, int lda )
|
||||
{
|
||||
int i, j;
|
||||
printf( "\n %s\n", desc );
|
||||
|
||||
for( i = 0; i < m; i++ )
|
||||
{
|
||||
for( j = 0; j < n; j++ ) printf( " %6.12f", a[i+j*lda] );
|
||||
printf( "\n" );
|
||||
}
|
||||
}
|
||||
|
||||
void MLPOD::print_matrix(const char *desc, int m, int n, int *a, int lda)
|
||||
{
|
||||
int i, j;
|
||||
printf( "\n %s\n", desc );
|
||||
|
||||
for( i = 0; i < m; i++ )
|
||||
{
|
||||
for( j = 0; j < n; j++ ) printf( " %d", a[i+j*lda] );
|
||||
printf( "\n" );
|
||||
}
|
||||
}
|
||||
|
||||
void MLPOD::podMatMul(double *c, double *a, double *b, int r1, int c1, int c2)
|
||||
{
|
||||
int i, j, k;
|
||||
@ -165,84 +128,6 @@ void MLPOD::podArrayFill(int* output, int start, int length)
|
||||
output[j] = start + j;
|
||||
}
|
||||
|
||||
double MLPOD::podArraySum(double *a, int n)
|
||||
{
|
||||
double e = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
e += a[i];
|
||||
return e;
|
||||
}
|
||||
|
||||
double MLPOD::podArrayMin(double *a, int n)
|
||||
{
|
||||
double b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]<b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
double MLPOD::podArrayMax(double *a, int n)
|
||||
{
|
||||
double b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]>b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
int MLPOD::podArrayMin(int *a, int n)
|
||||
{
|
||||
int b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]<b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
int MLPOD::podArrayMax(int *a, int n)
|
||||
{
|
||||
int b = a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
if (a[i]>b)
|
||||
b = a[i];
|
||||
return b;
|
||||
}
|
||||
|
||||
void MLPOD::podKron(double *C, double *A, double *B, double alpha, int M1, int M2)
|
||||
{
|
||||
int M = M1*M2;
|
||||
for (int idx=0; idx<M; idx++)
|
||||
{
|
||||
int ib = idx%M2;
|
||||
int ia = (idx-ib)/M2;
|
||||
C[idx] += alpha*A[ia]*B[ib];
|
||||
}
|
||||
}
|
||||
|
||||
void MLPOD::podCumsum(int* output, int* input, int length)
|
||||
{
|
||||
output[0] = 0;
|
||||
for (int j = 1; j < length; ++j)
|
||||
output[j] = input[j - 1] + output[j - 1];
|
||||
}
|
||||
|
||||
double MLPOD::podArrayNorm(double *a, int n)
|
||||
{
|
||||
double e = a[0]*a[0];
|
||||
for (int i=1; i<n; i++)
|
||||
e += a[i]*a[i];
|
||||
return sqrt(e);
|
||||
}
|
||||
|
||||
double MLPOD::podArrayErrorNorm(double *a, double *b, int n)
|
||||
{
|
||||
double e = (a[0]-b[0])*(a[0]-b[0]);
|
||||
for (int i=1; i<n; i++)
|
||||
e += (a[i]-b[i])*(a[i]-b[i]);
|
||||
return sqrt(e);
|
||||
}
|
||||
|
||||
void MLPOD::podArraySetValue(double *y, double a, int n)
|
||||
{
|
||||
for (int i=0; i<n; i++)
|
||||
@ -255,86 +140,6 @@ void MLPOD::podArrayCopy(double *y, double *x, int n)
|
||||
y[i] = x[i];
|
||||
}
|
||||
|
||||
void MLPOD::rotation_matrix(double *Rmat, double alpha, double beta, double gamma)
|
||||
{
|
||||
double ca = cos(alpha);
|
||||
double cb = cos(beta);
|
||||
double cg = cos(gamma);
|
||||
double sa = sin(alpha);
|
||||
double sb = sin(beta);
|
||||
double sg = sin(gamma);
|
||||
|
||||
Rmat[0] = ca*cg*cb - sa*sg;
|
||||
Rmat[3] = -ca*cb*sg - sa*cg;
|
||||
Rmat[6] = ca*sb;
|
||||
|
||||
Rmat[1] = sa*cg*cb + ca*sg;
|
||||
Rmat[4] = -sa*cb*sg + ca*cg;
|
||||
Rmat[7] = sa*sb;
|
||||
|
||||
Rmat[2] = -sb*cg;
|
||||
Rmat[5] = sb*sg;
|
||||
Rmat[8] = cb;
|
||||
}
|
||||
|
||||
void MLPOD::matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom)
|
||||
{
|
||||
double x1, x2, x3;
|
||||
for (int i=0; i < natom; i++) {
|
||||
x1 = x[0 + 3*i];
|
||||
x2 = x[1 + 3*i];
|
||||
x3 = x[2 + 3*i];
|
||||
xrot[0 + 3*i] = Rmat[0]*x1 + Rmat[3]*x2 + Rmat[6]*x3;
|
||||
xrot[1 + 3*i] = Rmat[1]*x1 + Rmat[4]*x2 + Rmat[7]*x3;
|
||||
xrot[2 + 3*i] = Rmat[2]*x1 + Rmat[5]*x2 + Rmat[8]*x3;
|
||||
}
|
||||
}
|
||||
|
||||
void MLPOD::matrix33_inverse(double *invA, double *A1, double *A2, double *A3)
|
||||
{
|
||||
double a11 = A1[0];
|
||||
double a21 = A1[1];
|
||||
double a31 = A1[2];
|
||||
double a12 = A2[0];
|
||||
double a22 = A2[1];
|
||||
double a32 = A2[2];
|
||||
double a13 = A3[0];
|
||||
double a23 = A3[1];
|
||||
double a33 = A3[2];
|
||||
double detA = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
|
||||
|
||||
invA[0] = (a22*a33 - a23*a32)/detA;
|
||||
invA[1] = (a23*a31 - a21*a33)/detA;
|
||||
invA[2] = (a21*a32 - a22*a31)/detA;
|
||||
invA[3] = (a13*a32 - a12*a33)/detA;
|
||||
invA[4] = (a11*a33 - a13*a31)/detA;
|
||||
invA[5] = (a12*a31 - a11*a32)/detA;
|
||||
invA[6] = (a12*a23 - a13*a22)/detA;
|
||||
invA[7] = (a13*a21 - a11*a23)/detA;
|
||||
invA[8] = (a11*a22 - a12*a21)/detA;
|
||||
}
|
||||
|
||||
void MLPOD::triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B, double *C)
|
||||
{
|
||||
double Anorm = sqrt(A[0]*A[0] + A[1]*A[1] + A[2]*A[2]);
|
||||
double Bnorm = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2]);
|
||||
double Cnorm = sqrt(C[0]*C[0] + C[1]*C[1] + C[2]*C[2]);
|
||||
|
||||
double Ahat[3];
|
||||
Ahat[0] = A[0]/Anorm; Ahat[1] = A[1]/Anorm; Ahat[2] = A[2]/Anorm;
|
||||
|
||||
double ax = Anorm;
|
||||
double bx = B[0]*Ahat[0] + B[1]*Ahat[1] + B[2]*Ahat[2]; //dot(B,Ahat);
|
||||
double by = sqrt(Bnorm*Bnorm - bx*bx); //sqrt(Bnorm^2 - bx^2);// #norm(cross(Ahat,B));
|
||||
double cx = C[0]*Ahat[0] + C[1]*Ahat[1] + C[2]*Ahat[2]; // dot(C,Ahat);
|
||||
double cy = (B[0]*C[0] + B[1]*C[1] + B[2]*C[2] - bx*cx)/by; // (dot(B, C) - bx*cx)/by;
|
||||
double cz = sqrt(Cnorm*Cnorm - cx*cx - cy*cy); // sqrt(Cnorm^2 - cx^2 - cy^2);
|
||||
|
||||
a[0] = ax; a[1] = 0.0; a[2] = 0.0;
|
||||
b[0] = bx; b[1] = by; b[2] = 0.0;
|
||||
c[0] = cx; c[1] = cy; c[2] = cz;
|
||||
}
|
||||
|
||||
void podsnapshots(double *rbf, double *xij, double *besselparams, double rin, double rcut,
|
||||
int besseldegree, int inversedegree, int nbesselpars, int N)
|
||||
{
|
||||
|
||||
@ -253,27 +253,10 @@ class MLPOD : protected Pointers {
|
||||
|
||||
// functions for collecting/collating arrays
|
||||
|
||||
void print_matrix(const char *desc, int m, int n, int *a, int lda);
|
||||
void print_matrix(const char *desc, int m, int n, double *a, int lda);
|
||||
void print_matrix(const char *desc, int m, int n, double **a, int lda);
|
||||
void podMatMul(double *c, double *a, double *b, int r1, int c1, int c2);
|
||||
void podCumsum(int *output, int *input, int length);
|
||||
double podArrayNorm(double *a, int n);
|
||||
double podArrayErrorNorm(double *a, double *b, int n);
|
||||
void podArraySetValue(double *y, double a, int n);
|
||||
void podArrayCopy(double *y, double *x, int n);
|
||||
void podArrayFill(int *output, int start, int length);
|
||||
double podArrayMin(double *a, int n);
|
||||
double podArrayMax(double *a, int n);
|
||||
double podArraySum(double *a, int n);
|
||||
int podArrayMin(int *a, int n);
|
||||
int podArrayMax(int *a, int n);
|
||||
void podKron(double *C, double *A, double *B, double alpha, int M1, int M2);
|
||||
void rotation_matrix(double *Rmat, double alpha, double beta, double gamma);
|
||||
void triclinic_lattice_conversion(double *a, double *b, double *c, double *A, double *B,
|
||||
double *C);
|
||||
void matrix33_multiplication(double *xrot, double *Rmat, double *x, int natom);
|
||||
void matrix33_inverse(double *invA, double *A1, double *A2, double *A3);
|
||||
|
||||
// functions for calculating energy and force descriptors
|
||||
|
||||
|
||||
Reference in New Issue
Block a user