diff --git a/doc/src/compute_mliap.rst b/doc/src/compute_mliap.rst index 91253b29c2..3eb2a8832a 100644 --- a/doc/src/compute_mliap.rst +++ b/doc/src/compute_mliap.rst @@ -20,7 +20,7 @@ Syntax *model* values = style Nelems Nparams style = *linear* or *quadratic* Nelems = number of elements - Nparams = number of parameters per element + Nparams = number of parameters per element *descriptor* values = style filename style = *sna* filename = name of file containing descriptor definitions @@ -36,25 +36,25 @@ Description """"""""""" Compute style *mliap* provides a general interface to the gradient -of machine-learning interatomic potentials w.r.t. model parameters. +of machine-learning interatomic potentials w.r.t. model parameters. It is used primarily for calculating the gradient of energy, force, and stress components w.r.t. model parameters, which is useful when training :doc:`mliap pair_style ` models to match target data. -It provides separate +It provides separate definitions of the interatomic potential functional form (*model*) and the geometric quantities that characterize the atomic positions -(*descriptor*). By defining *model* and *descriptor* separately, +(*descriptor*). By defining *model* and *descriptor* separately, it is possible to use many different models with a given descriptor, -or many different descriptors with a given model. Currently, the +or many different descriptors with a given model. Currently, the compute supports just two models, *linear* and *quadratic*, -and one descriptor, *sna*, the SNAP descriptor used by +and one descriptor, *sna*, the SNAP descriptor used by :doc:`pair_style snap `, including the linear, quadratic, and chem variants. Work is currently underway to extend the interface to handle neural network energy models, and it is also straightforward to add new descriptor styles. The compute *mliap* command must be followed by two keywords -*model* and *descriptor* in either order. +*model* and *descriptor* in either order. The *model* keyword is followed by a model style, currently limited to either *linear* or *quadratic*. In both cases, @@ -62,25 +62,25 @@ this is followed by two arguments. *nelems* is the number of elements. It must be equal to the number of LAMMPS atom types. *nparams* is the number of parameters per element for this model i.e. the number of parameter gradients for each element. Note these definitions -are identical to those of *nelems* and *nparams* in the +are identical to those of *nelems* and *nparams* in the :doc:`pair_style mliap ` model file. - + The *descriptor* keyword is followed by a descriptor style, and additional arguments. -Currently the only descriptor style is *sna*, indicating the bispectrum component -descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of +Currently the only descriptor style is *sna*, indicating the bispectrum component +descriptors used by the Spectral Neighbor Analysis Potential (SNAP) potentials of :doc:`pair_style snap `. -The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped -out by the LAMMPS command parser. A single additional argument specifies the descriptor filename -containing the parameters and setting used by the SNAP descriptor. +The \'p\' in SNAP is dropped, because keywords that match pair_styles are silently stripped +out by the LAMMPS command parser. A single additional argument specifies the descriptor filename +containing the parameters and setting used by the SNAP descriptor. The descriptor filename usually ends in the *.mliap.descriptor* extension. -The format of this file is identical to the descriptor file in the +The format of this file is identical to the descriptor file in the :doc:`pair_style mliap `, and is described in detail -there. +there. .. note:: The number of LAMMPS atom types (and the value of *nelems* in the model) - must match the value of *nelems* in the descriptor file. + must match the value of *nelems* in the descriptor file. Compute *mliap* calculates a global array containing gradient information. The number of columns in the array is :math:`nelems \times nparams + 1`. @@ -89,9 +89,9 @@ each parameter and each element. The last six rows of the array contain the corresponding derivatives of the virial stress tensor, listed in Voigt notation: *pxx*, *pyy*, *pzz*, *pyz*, *pxz*, *pxy*. In between the energy and stress rows are -the 3\*\ *N* rows containing the derivatives of the force components. -See section below on output for a detailed description of how -rows and columns are ordered. +the 3\*\ *N* rows containing the derivatives of the force components. +See section below on output for a detailed description of how +rows and columns are ordered. The element in the last column of each row contains the potential energy, force, or stress, according to the row. @@ -109,7 +109,7 @@ command: See section below on output for a detailed explanation of the data layout in the global array. -Atoms not in the group do not contribute to this compute. +Atoms not in the group do not contribute to this compute. Neighbor atoms not in the group do not contribute to this compute. The neighbor list needed to compute this quantity is constructed each time the calculation is performed (i.e. each time a snapshot of atoms @@ -120,7 +120,7 @@ too frequently. If the user-specified reference potentials includes bonded and non-bonded pairwise interactions, then the settings of - :doc:`special_bonds ` command can remove pairwise + :doc:`special_bonds ` command can remove pairwise interactions between atoms in the same bond, angle, or dihedral. This is the default setting for the :doc:`special_bonds ` command, and means those pairwise interactions do not appear in the @@ -128,7 +128,7 @@ too frequently. those pairs will not be included in the calculation. The :doc:`rerun ` command is not an option here, since the reference potential is required for the last column of the global array. A work-around is to prevent - pairwise interactions from being removed by explicitly adding a + pairwise interactions from being removed by explicitly adding a *tiny* positive value for every pairwise interaction that would otherwise be set to zero in the :doc:`special_bonds ` command. @@ -139,7 +139,7 @@ too frequently. Compute *mliap* evaluates a global array. The columns are arranged into *nelems* blocks, listed in order of element *I*\ . Each block -contains one column for each of the *nparams* model parameters. +contains one column for each of the *nparams* model parameters. A final column contains the corresponding energy, force component on an atom, or virial stress component. The rows of the array appear in the following order: diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index cf9339a30d..4c8c7e1383 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -392,7 +392,7 @@ of :math:`K N_{elem}^3` columns. These values can be accessed by any command that uses per-atom values from a compute as input. See the :doc:`Howto output ` doc page for an overview of LAMMPS output options. To see how this command -can be used within a Python workflow to train SNAP potentials, +can be used within a Python workflow to train SNAP potentials, see the examples in `FitSNAP `_. Restrictions diff --git a/doc/src/pair_lj_cut_coul.rst b/doc/src/pair_lj_cut_coul.rst index fda3fd29c5..38553a205f 100644 --- a/doc/src/pair_lj_cut_coul.rst +++ b/doc/src/pair_lj_cut_coul.rst @@ -7,20 +7,20 @@ .. index:: pair_style lj/cut/coul/debye/kk .. index:: pair_style lj/cut/coul/debye/omp .. index:: pair_style lj/cut/coul/dsf -.. index:: pair_style lj/cut/coul/dsf/gpu -.. index:: pair_style lj/cut/coul/dsf/kk -.. index:: pair_style lj/cut/coul/dsf/omp -.. index:: pair_style lj/cut/coul/long -.. index:: pair_style lj/cut/coul/long/gpu -.. index:: pair_style lj/cut/coul/long/kk -.. index:: pair_style lj/cut/coul/long/intel -.. index:: pair_style lj/cut/coul/long/opt -.. index:: pair_style lj/cut/coul/long/omp -.. index:: pair_style lj/cut/coul/msm -.. index:: pair_style lj/cut/coul/msm/gpu -.. index:: pair_style lj/cut/coul/msm/omp -.. index:: pair_style lj/cut/coul/wolf -.. index:: pair_style lj/cut/coul/wolf/omp +.. index:: pair_style lj/cut/coul/dsf/gpu +.. index:: pair_style lj/cut/coul/dsf/kk +.. index:: pair_style lj/cut/coul/dsf/omp +.. index:: pair_style lj/cut/coul/long +.. index:: pair_style lj/cut/coul/long/gpu +.. index:: pair_style lj/cut/coul/long/kk +.. index:: pair_style lj/cut/coul/long/intel +.. index:: pair_style lj/cut/coul/long/opt +.. index:: pair_style lj/cut/coul/long/omp +.. index:: pair_style lj/cut/coul/msm +.. index:: pair_style lj/cut/coul/msm/gpu +.. index:: pair_style lj/cut/coul/msm/omp +.. index:: pair_style lj/cut/coul/wolf +.. index:: pair_style lj/cut/coul/wolf/omp pair_style lj/cut/coul/cut command ================================== @@ -234,7 +234,7 @@ LJ and Coulombic cutoffs specified in the pair_style command are used. If only one cutoff is specified, it is used as the cutoff for both LJ and Coulombic interactions for this type pair. If both coefficients are specified, they are used as the LJ and Coulombic cutoffs for this -type pair. +type pair. For *lj/cut/coul/long* and *lj/cut/coul/msm* only the LJ cutoff can be specified since a Coulombic cutoff cannot be specified for an individual I,J diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index d96e35a7f4..c93dc4981f 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -34,8 +34,8 @@ Examples Description """"""""""" -Pair style *mliap* provides a general interface to families of -machine-learning interatomic potentials. It allows separate +Pair style *mliap* provides a general interface to families of +machine-learning interatomic potentials. It allows separate definitions of the interatomic potential functional form (*model*) and the geometric quantities that characterize the atomic positions (*descriptor*). By defining *model* and *descriptor* separately, @@ -61,8 +61,8 @@ where N is the number of LAMMPS atom types. The *model* keyword is followed by a model style, currently limited to either *linear* or *quadratic*. In both cases, -this is followed by a single argument specifying the model filename containing the -parameters for a set of elements. +this is followed by a single argument specifying the model filename containing the +parameters for a set of elements. The model filename usually ends in the *.mliap.model* extension. It may contain parameters for many elements. The only requirement is that it contain at least those element names appearing in the diff --git a/src/BODY/body_nparticle.cpp b/src/BODY/body_nparticle.cpp index a3ffc9a5b5..21fb2a72ed 100644 --- a/src/BODY/body_nparticle.cpp +++ b/src/BODY/body_nparticle.cpp @@ -213,7 +213,7 @@ int BodyNparticle::pack_data_body(tagint atomID, int ibonus, double *buf) int nsub = ivalue[0]; if (buf) { - + // ID ninteger ndouble m = 0; @@ -222,7 +222,7 @@ int BodyNparticle::pack_data_body(tagint atomID, int ibonus, double *buf) buf[m++] = ubuf(6 + 3*nsub).d; // single integer Nsub - + buf[m++] = ubuf(nsub).d; // 6 moments of inertia @@ -239,7 +239,7 @@ int BodyNparticle::pack_data_body(tagint atomID, int ibonus, double *buf) buf[m++] = ispace[1][2]; // 3*Nsub particle coords = displacement from COM in box frame - + for (int i = 0; i < nsub; i++) { MathExtra::matvec(p,&dvalue[3*i],values); buf[m++] = values[0]; @@ -248,7 +248,7 @@ int BodyNparticle::pack_data_body(tagint atomID, int ibonus, double *buf) } } else m = 3 + 1 + 6 + 3*nsub; - + return m; } @@ -272,10 +272,10 @@ int BodyNparticle::write_data_body(FILE *fp, double *buf) fmt::print(fp,"{} {} {} {} {} {}\n", inertia[0],inertia[1],inertia[2],inertia[3],inertia[4],inertia[5]); m += 6; - + for (int i = 0; i < nsub; i++) fmt::print(fp,"{} {} {}\n",buf[m++],buf[m++],buf[m++]); - + return m; } diff --git a/src/MLIAP/README b/src/MLIAP/README index 77e38d3c62..18ce3297e6 100644 --- a/src/MLIAP/README +++ b/src/MLIAP/README @@ -1,16 +1,16 @@ -This package provides a general interface to families of +This package provides a general interface to families of machine-learning interatomic potentials (MLIAPs). This interface consists -of a mliap pair style and a mliap compute. +of a mliap pair style and a mliap compute. -The mliap pair style is used when running simulations with energies and -forces calculated by an MLIAP. The interface allows separate +The mliap pair style is used when running simulations with energies and +forces calculated by an MLIAP. The interface allows separate definitions of the interatomic potential functional form (*model*) and the geometric quantities that characterize the atomic positions -(*descriptor*). By defining *model* and *descriptor* separately, +(*descriptor*). By defining *model* and *descriptor* separately, it is possible to use many different models with a given descriptor, -or many different descriptors with a given model. Currently, the pair_style +or many different descriptors with a given model. Currently, the pair_style supports just two models, *linear* and *quadratic*, -and one descriptor, *sna*, the SNAP descriptor, including the +and one descriptor, *sna*, the SNAP descriptor, including the linear, quadratic, and chem variants. Work is currently underway to extend the interface to handle neural network energy models, and it is also straightforward to add new descriptor styles. @@ -22,8 +22,8 @@ Any *model or *descriptor* that has been implemented for the *mliap* pair style can also be accessed by the *mliap* compute. In addition to the energy, force, and stress gradients, w.r.t. each *model* parameter, the compute also calculates the energy, -force, and stress contributions from a user-specified -reference potential. +force, and stress contributions from a user-specified +reference potential. To see how this command can be used within a Python workflow to train machine-learning interatomic diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp index 655e70c6f9..77bae76d7b 100644 --- a/src/MLIAP/compute_mliap.cpp +++ b/src/MLIAP/compute_mliap.cpp @@ -35,7 +35,7 @@ enum{SCALAR,VECTOR,ARRAY}; ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), list(NULL), mliap(NULL), - gradforce(NULL), mliapall(NULL), map(NULL), + gradforce(NULL), mliapall(NULL), map(NULL), descriptors(NULL), gamma_row_index(NULL), gamma_col_index(NULL), gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL) { @@ -263,15 +263,15 @@ void ComputeMLIAP::compute_array() // i.e. gamma = d2E/dsigma_l.dB_k // sigma_l is a parameter and B_k is a descriptor of atom i // for SNAP, this is a sparse natoms*nparams*ndescriptors matrix, - // but in general it could be fully dense. - - model->param_gradient(map, list, descriptors, gamma_row_index, + // but in general it could be fully dense. + + model->param_gradient(map, list, descriptors, gamma_row_index, gamma_col_index, gamma, egradient); // calculate descriptor gradient contributions to parameter gradients - descriptor->compute_gradients(map, list, gamma_nnz, gamma_row_index, + descriptor->compute_gradients(map, list, gamma_nnz, gamma_row_index, gamma_col_index, gamma, gradforce, yoffset, zoffset); diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h index 74b5fc1f93..d25b5e66b2 100644 --- a/src/MLIAP/compute_mliap.h +++ b/src/MLIAP/compute_mliap.h @@ -44,13 +44,13 @@ class ComputeMLIAP : public Compute { int nelements; double** descriptors; // descriptors for all atoms in list - int ndescriptors; // number of descriptors + int ndescriptors; // number of descriptors int gamma_max; // number of atoms allocated for beta, descriptors int nparams; // number of model paramters per element int gamma_nnz; // number of non-zero entries in gamma double** gamma; // gamma element - int** gamma_row_index; // row (parameter) index - int** gamma_col_index; // column (descriptor) index + int** gamma_row_index; // row (parameter) index + int** gamma_col_index; // column (descriptor) index double* egradient; // energy gradient w.r.t. parameters class MLIAPModel* model; diff --git a/src/MLIAP/mliap_descriptor.h b/src/MLIAP/mliap_descriptor.h index 072a5b4729..c6f06c52db 100644 --- a/src/MLIAP/mliap_descriptor.h +++ b/src/MLIAP/mliap_descriptor.h @@ -24,9 +24,9 @@ public: ~MLIAPDescriptor(); virtual void compute_descriptors(int*, class NeighList*, double**)=0; virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int)=0; - virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**, + virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**, double**, int, int)=0; - virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**, + virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**, double**, int, int)=0; virtual void init()=0; virtual double memory_usage()=0; @@ -34,7 +34,7 @@ public: int ndescriptors; // number of descriptors int nelements; // # of unique elements char **elements; // names of unique elements - double **cutsq; // nelem x nelem rcutsq values + double **cutsq; // nelem x nelem rcutsq values double cutmax; // maximum cutoff needed protected: diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index c2879e942a..f11c80e04f 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -259,8 +259,8 @@ void MLIAPDescriptorSNAP::compute_forces(PairMLIAP* pairmliap, NeighList* list, compute force gradient for each atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, - int gamma_nnz, int **gamma_row_index, +void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, + int gamma_nnz, int **gamma_row_index, int **gamma_col_index, double **gamma, double **gradforce, int yoffset, int zoffset) { @@ -345,9 +345,9 @@ void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, snaptr->rcutij[jj],jj, 0); snaptr->compute_dbidrj(); - + // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj - + for (int inz = 0; inz < gamma_nnz; inz++) { const int l = gamma_row_index[ii][inz]; const int k = gamma_col_index[ii][inz]; @@ -358,7 +358,7 @@ void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, gradforce[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1]; gradforce[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2]; } - + } } @@ -368,8 +368,8 @@ void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, compute descriptor gradients for each neighbor atom ---------------------------------------------------------------------- */ -void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list, - int gamma_nnz, int **gamma_row_index, +void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list, + int gamma_nnz, int **gamma_row_index, int **gamma_col_index, double **gamma, double **graddesc, int yoffset, int zoffset) { @@ -464,7 +464,7 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list graddesc[j][k] = -snaptr->dblist[k][0]; graddesc[j][k] = -snaptr->dblist[k][1]; graddesc[j][k] = -snaptr->dblist[k][2]; - } + } } } diff --git a/src/MLIAP/mliap_descriptor_snap.h b/src/MLIAP/mliap_descriptor_snap.h index a929e01149..eda5a8f111 100644 --- a/src/MLIAP/mliap_descriptor_snap.h +++ b/src/MLIAP/mliap_descriptor_snap.h @@ -24,9 +24,9 @@ public: ~MLIAPDescriptorSNAP(); virtual void compute_descriptors(int*, class NeighList*, double**); virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int); - virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**, + virtual void compute_gradients(int*, class NeighList*, int, int**, int**, double**, double**, int, int); - virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**, + virtual void compute_descriptor_gradients(int*, class NeighList*, int, int**, int**, double**, double**, int, int); virtual void init(); virtual double memory_usage(); diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp index b8d06155ba..b9358f9f75 100644 --- a/src/MLIAP/mliap_model_linear.cpp +++ b/src/MLIAP/mliap_model_linear.cpp @@ -36,7 +36,7 @@ MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) : /* ---------------------------------------------------------------------- */ -MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in) : +MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in) : MLIAPModel(lmp, nelements_in, nparams_in) { ndescriptors = nparams - 1; @@ -47,7 +47,7 @@ MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in MLIAPModelLinear::~MLIAPModelLinear(){} /* ---------------------------------------------------------------------- - Calculate model gradients w.r.t descriptors + Calculate model gradients w.r.t descriptors for each atom beta_i = dE(B_i)/dB_i ---------------------------------------------------------------------- */ @@ -87,22 +87,22 @@ void MLIAPModelLinear::gradient(PairMLIAP* pairmliap, NeighList* list, double ** /* ---------------------------------------------------------------------- Calculate model double gradients w.r.t descriptors and parameters - for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l, - where sigma_l is a parameter, B_k a descriptor, + for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l, + where sigma_l is a parameter, B_k a descriptor, and atom subscript i is omitted gamma is in CSR format: nnz = number of non-zero values - gamma_row_index[inz] = l indices, 0 <= l < nparams + gamma_row_index[inz] = l indices, 0 <= l < nparams gamma_col_indexiinz] = k indices, 0 <= k < ndescriptors gamma[i][inz] = non-zero values, 0 <= inz < nnz egradient is derivative of energy w.r.t. parameters ---------------------------------------------------------------------- */ -void MLIAPModelLinear::param_gradient(int *map, NeighList* list, - double **descriptors, - int **gamma_row_index, int **gamma_col_index, +void MLIAPModelLinear::param_gradient(int *map, NeighList* list, + double **descriptors, + int **gamma_row_index, int **gamma_col_index, double **gamma, double *egradient) { int i; @@ -112,7 +112,7 @@ void MLIAPModelLinear::param_gradient(int *map, NeighList* list, for (int l = 0; l < nelements*nparams; l++) egradient[l] = 0.0; - + for (int ii = 0; ii < list->inum; ii++) { i = list->ilist[ii]; @@ -128,12 +128,12 @@ void MLIAPModelLinear::param_gradient(int *map, NeighList* list, } // gradient of energy of atom I w.r.t. parameters - + l = elemoffset; egradient[l++] += 1.0; for (int icoeff = 0; icoeff < ndescriptors; icoeff++) egradient[l++] += descriptors[ii][icoeff]; - + } } diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp index f5aae1d89d..8e8f747628 100644 --- a/src/MLIAP/mliap_model_quadratic.cpp +++ b/src/MLIAP/mliap_model_quadratic.cpp @@ -37,7 +37,7 @@ MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) : /* ---------------------------------------------------------------------- */ -MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, int nelements_in, int nparams_in) : +MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, int nelements_in, int nparams_in) : MLIAPModel(lmp, nelements_in, nparams_in) { nonlinearflag = 1; @@ -113,22 +113,22 @@ void MLIAPModelQuadratic::gradient(PairMLIAP* pairmliap, NeighList* list, double /* ---------------------------------------------------------------------- Calculate model double gradients w.r.t descriptors and parameters - for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l, - where sigma_l is a parameter, B_k a descriptor, + for each atom energy gamma_lk = d2E(B)/dB_k/dsigma_l, + where sigma_l is a parameter, B_k a descriptor, and atom subscript i is omitted gamma is in CSR format: nnz = number of non-zero values - gamma_row_index[inz] = l indices, 0 <= l < nparams + gamma_row_index[inz] = l indices, 0 <= l < nparams gamma_col_indexiinz] = k indices, 0 <= k < ndescriptors gamma[i][inz] = non-zero values, 0 <= inz < nnz egradient is derivative of energy w.r.t. parameters ---------------------------------------------------------------------- */ -void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list, - double **descriptors, - int **gamma_row_index, int **gamma_col_index, +void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list, + double **descriptors, + int **gamma_row_index, int **gamma_col_index, double **gamma, double *egradient) { int i; @@ -138,7 +138,7 @@ void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list, for (int l = 0; l < nelements*nparams; l++) egradient[l] = 0.0; - + for (int ii = 0; ii < list->inum; ii++) { i = list->ilist[ii]; @@ -179,14 +179,14 @@ void MLIAPModelQuadratic::param_gradient(int *map, NeighList* list, } // gradient of energy of atom I w.r.t. parameters - + l = elemoffset; egradient[l++] += 1.0; for (int icoeff = 0; icoeff < ndescriptors; icoeff++) egradient[l++] += descriptors[ii][icoeff]; - + // quadratic contributions - + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { double bveci = descriptors[ii][icoeff]; egradient[l++] += 0.5*bveci*bveci; diff --git a/src/USER-MESONT/pair_mesocnt.h b/src/USER-MESONT/pair_mesocnt.h index 770aa21e36..2e13fe9243 100644 --- a/src/USER-MESONT/pair_mesocnt.h +++ b/src/USER-MESONT/pair_mesocnt.h @@ -34,7 +34,7 @@ class PairMesoCNT : public Pair { double r_ang,rsq_ang,d_ang,rc_ang,rcsq_ang,cutoff_ang,cutoffsq_ang; double sig,sig_ang,comega,ctheta; double hstart_uinf,hstart_gamma, - hstart_phi,psistart_phi,hstart_usemi,xistart_usemi; + hstart_phi,psistart_phi,hstart_usemi,xistart_usemi; double delh_uinf,delh_gamma,delh_phi,delpsi_phi,delh_usemi,delxi_usemi; double p1[3],p2[3],p[3],m[3]; @@ -55,7 +55,7 @@ class PairMesoCNT : public Pair { void read_file(); void read_data(FILE *, double *, double &, double &, int); void read_data(FILE *, double **, double &, double &, - double &, double &, int); + double &, double &, int); void spline_coeff(double *, double **, double, int); void spline_coeff(double **, double ****, double, double, int); @@ -63,17 +63,17 @@ class PairMesoCNT : public Pair { double spline(double, double, double, double **, int); double dspline(double, double, double, double **, int); double spline(double, double, double, double, double, double, - double ****, int); + double ****, int); double dxspline(double, double, double, double, double, double, - double ****, int); + double ****, int); double dyspline(double, double, double, double, double, double, - double ****, int); + double ****, int); void geometry(const double *, const double *, const double *, - const double *, const double *, + const double *, const double *, double *, double *, double *, double **); void weight(const double *, const double *, const double *, - const double *, double &, double *, double *, + const double *, double &, double *, double *, double *, double *); void finf(const double *, double &, double **); @@ -97,7 +97,7 @@ class PairMesoCNT : public Pair { inline double s5(double x) { double x2 = x * x; return heaviside(-x) - + heaviside(x)*heaviside(1-x)*(1 - x2*x*(6*x2 - 15*x + 10)); + + heaviside(x)*heaviside(1-x)*(1 - x2*x*(6*x2 - 15*x + 10)); } inline double ds5(double x) { diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 04522882f9..c944cf3eda 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -587,7 +587,7 @@ void AtomVecBody::pack_data_pre(int ilocal) int AtomVecBody::pack_data_bonus(double *buf, int /*flag*/) { int i,j; - + tagint *tag = atom->tag; int nlocal = atom->nlocal; @@ -596,7 +596,7 @@ int AtomVecBody::pack_data_bonus(double *buf, int /*flag*/) if (body[i] < 0) continue; m += bptr->pack_data_body(tag[i],body[i],buf); } - + return m; } diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index ac7737bcfc..a1b8079192 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -490,7 +490,7 @@ void AtomVecEllipsoid::pack_data_post(int ilocal) int AtomVecEllipsoid::pack_data_bonus(double *buf, int /*flag*/) { int i,j; - + tagint *tag = atom->tag; int nlocal = atom->nlocal; diff --git a/src/atom_vec_ellipsoid.h b/src/atom_vec_ellipsoid.h index a7961f97f8..041a20affd 100644 --- a/src/atom_vec_ellipsoid.h +++ b/src/atom_vec_ellipsoid.h @@ -58,7 +58,7 @@ class AtomVecEllipsoid : public AtomVec { int pack_data_bonus(double *, int); void write_data_bonus(FILE *, int, double *, int); - + // unique to AtomVecEllipsoid void set_shape(int, double, double, double); diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index f024ed8a00..db593dcee0 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -469,7 +469,7 @@ int AtomVecHybrid::pack_data_bonus(double *buf, int flag) if (flag == LINE && strcmp(keywords[k],"line") != 0) continue; if (flag == TRIANGLE && strcmp(keywords[k],"tri") != 0) continue; if (flag == BODY && strcmp(keywords[k],"body") != 0) continue; - + return styles[k]->pack_data_bonus(buf,flag); } return 0; diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp index f6e42933fc..26779811e0 100644 --- a/src/atom_vec_line.cpp +++ b/src/atom_vec_line.cpp @@ -465,7 +465,7 @@ int AtomVecLine::pack_data_bonus(double *buf, int /*flag*/) int i,j; double length,theta; double xc,yc,x1,x2,y1,y2; - + double **x = atom->x; tagint *tag = atom->tag; int nlocal = atom->nlocal; diff --git a/src/body.h b/src/body.h index e843c21b40..1c742d3e65 100644 --- a/src/body.h +++ b/src/body.h @@ -45,7 +45,7 @@ class Body : protected Pointers { virtual void data_body(int, int, int, int*, double *) = 0; virtual int pack_data_body(tagint, int, double *) = 0; virtual int write_data_body(FILE *, double *) = 0; - + virtual int noutrow(int) = 0; virtual int noutcol() = 0; virtual void output(int, int, double *) = 0; diff --git a/src/write_data.cpp b/src/write_data.cpp index 60d3781ff5..0804387577 100644 --- a/src/write_data.cpp +++ b/src/write_data.cpp @@ -198,14 +198,14 @@ void WriteData::write(const std::string &file) // molecular topology info if defined // do not write molecular topology for atom_style template - + if (atom->molecular == 1) { if (atom->nbonds && nbonds) bonds(); if (atom->nangles && nangles) angles(); if (atom->ndihedrals) dihedrals(); if (atom->nimpropers) impropers(); } - + // bonus info if defined if (natoms && atom->ellipsoid_flag) bonus(ELLIPSOID); @@ -214,7 +214,7 @@ void WriteData::write(const std::string &file) if (natoms && atom->body_flag) bonus(BODY); // extra sections managed by fixes - + if (fixflag) for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->wd_section) @@ -261,7 +261,7 @@ void WriteData::header() if (atom->body_flag) fmt::print(fp,"{} bodies\n",atom->nbodies); // fix info - + if (fixflag) for (int i = 0; i < modify->nfix; i++) if (modify->fix[i]->wd_header) @@ -671,7 +671,7 @@ void WriteData::bonus(int flag) // pack my bonus data into buf atom->avec->pack_data_bonus(buf,flag); - + // write one chunk of info per proc to file // proc 0 pings each proc, receives its chunk, writes to file // all other procs wait for ping, send their chunk to proc 0 @@ -686,7 +686,7 @@ void WriteData::bonus(int flag) if (flag == LINE) fprintf(fp,"\nLines\n\n"); if (flag == TRIANGLE) fprintf(fp,"\nTriangles\n\n"); if (flag == BODY) fprintf(fp,"\nBodies\n\n"); - + for (int iproc = 0; iproc < nprocs; iproc++) { if (iproc) { MPI_Irecv(buf,maxvalues,MPI_DOUBLE,iproc,0,world,&request);