diff --git a/examples/mliap/in.mliap.quadratic.compute b/examples/mliap/in.mliap.quadratic.compute index 6b346392fa..deaab169f1 100644 --- a/examples/mliap/in.mliap.quadratic.compute +++ b/examples/mliap/in.mliap.quadratic.compute @@ -69,7 +69,7 @@ variable db_2_100 equal c_db[2][100] # set up compute snap generating global array -compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 +compute snap all mliap descriptor sna compute.mliap.descriptor model quadratic 2 21 gradgradflag 0 fix snap all ave/time 1 1 1 c_snap[*] file compute.quadratic.dat mode vector thermo 100 diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute index 4b88a3e52d..1f638a61bb 100644 --- a/examples/mliap/in.mliap.snap.compute +++ b/examples/mliap/in.mliap.snap.compute @@ -69,7 +69,7 @@ variable db_2_25 equal c_db[2][25] # set up compute snap generating global array -compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 +compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 gradgradflag 1 fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector thermo 100 diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp index 655e70c6f9..7588248a6f 100644 --- a/src/MLIAP/compute_mliap.cpp +++ b/src/MLIAP/compute_mliap.cpp @@ -37,7 +37,9 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), list(NULL), mliap(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) + gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL), + iatomdesc(NULL), ielemdesc(NULL), numneighdesc(NULL), + jatomdesc(NULL), jelemdesc(NULL), graddesc(NULL) { array_flag = 1; extarray = 0; @@ -45,6 +47,10 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : if (narg < 4) error->all(FLERR,"Illegal compute mliap command"); + // default values + + gradgradflag = 1; + // set flags for required keywords int modelflag = 0; @@ -79,6 +85,12 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : iarg += 3; } else error->all(FLERR,"Illegal compute mliap command"); descriptorflag = 1; + } else if (strcmp(arg[iarg],"gradgradflag") == 0) { + if (iarg+1 > narg) error->all(FLERR,"Illegal compute mliap command"); + gradgradflag = atoi(arg[iarg+1]); + if (gradgradflag != 0 && gradgradflag != 1) + error->all(FLERR,"Illegal compute mliap command"); + iarg += 2; } else error->all(FLERR,"Illegal compute mliap command"); } @@ -102,7 +114,8 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : size_gradforce = ndims_force*nparams*nelements; nmax = 0; - gamma_max = 0; + numlistdesc_max = 0; + nneighdesc_max = 0; // create a minimal map, placeholder for more general map @@ -129,6 +142,13 @@ ComputeMLIAP::~ComputeMLIAP() memory->destroy(gamma); memory->destroy(egradient); + memory->destroy(iatomdesc); + memory->destroy(ielemdesc); + memory->destroy(numneighdesc); + memory->destroy(jatomdesc); + memory->destroy(jelemdesc); + memory->destroy(graddesc); + delete model; delete descriptor; } @@ -245,35 +265,139 @@ void ComputeMLIAP::compute_array() neighbor->build_one(list); - if (gamma_max < list->inum) { - memory->grow(descriptors,list->inum,ndescriptors,"ComputeMLIAP:descriptors"); - memory->grow(gamma_row_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_row_index"); - memory->grow(gamma_col_index,list->inum,gamma_nnz,"ComputeMLIAP:gamma_col_index"); - memory->grow(gamma,list->inum,gamma_nnz,"ComputeMLIAP:gamma"); - gamma_max = list->inum; + const int numlistdesc = list->inum; + if (numlistdesc_max < numlistdesc) { + memory->grow(descriptors,numlistdesc,ndescriptors,"ComputeMLIAP:descriptors"); + memory->grow(gamma_row_index,numlistdesc,gamma_nnz,"ComputeMLIAP:gamma_row_index"); + memory->grow(gamma_col_index,numlistdesc,gamma_nnz,"ComputeMLIAP:gamma_col_index"); + memory->grow(gamma,numlistdesc,gamma_nnz,"ComputeMLIAP:gamma"); + memory->grow(iatomdesc,numlistdesc,"ComputeMLIAP:iatomdesc"); + memory->grow(ielemdesc,numlistdesc,"ComputeMLIAP:ielemdesc"); + memory->grow(numneighdesc,numlistdesc,"ComputeMLIAP:numneighdesc"); + numlistdesc_max = numlistdesc; } - // compute descriptors + if (gradgradflag) { - descriptor->compute_descriptors(map, list, descriptors); + // compute descriptors + + descriptor->compute_descriptors(map, list, descriptors); + + // calculate descriptor contributions to parameter gradients + // and gamma = double gradient w.r.t. parameters and descriptors + + // 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, + gamma_col_index, gamma, egradient); + + + // calculate descriptor gradient contributions to parameter gradients + + descriptor->compute_gradients(map, list, gamma_nnz, gamma_row_index, + gamma_col_index, gamma, gradforce, + yoffset, zoffset); + + } else { - // calculate descriptor contributions to parameter gradients - // and gamma = double gradient w.r.t. parameters and descriptors - - // 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, - gamma_col_index, gamma, egradient); - - - // calculate descriptor gradient contributions to parameter gradients - - descriptor->compute_gradients(map, list, gamma_nnz, gamma_row_index, - gamma_col_index, gamma, gradforce, - yoffset, zoffset); + double **x = atom->x; + int *type = atom->type; + + int *numneigh = list->numneigh; + int **firstneigh = list->firstneigh; + + int nneighdesc = 0; + for (int ii = 0; ii < list->inum; ii++) { + const int i = list->ilist[ii]; + + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const int ielem = map[itype]; + + int *jlist = firstneigh[i]; + const int jnum = numneigh[i]; + + int ninside = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + j &= NEIGHMASK; + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type[j]; + const int jelem = map[jtype]; + + if (rsq < descriptor->cutsq[ielem][jelem]) ninside++; + + } + iatomdesc[ii] = i; + ielemdesc[ii] = ielem; + numneighdesc[ii] = ninside; + nneighdesc += ninside; + } + + if (nneighdesc_max < nneighdesc) { + memory->grow(jatomdesc,nneighdesc,"ComputeMLIAP:jatomdesc"); + memory->grow(jelemdesc,nneighdesc,"ComputeMLIAP:jelemdesc"); + memory->grow(graddesc,nneighdesc,ndescriptors,3,"ComputeMLIAP:graddesc"); + nneighdesc_max = nneighdesc; + } + + int ij = 0; + for (int ii = 0; ii < list->inum; ii++) { + const int i = list->ilist[ii]; + + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + const int itype = type[i]; + const int ielem = map[itype]; + + int *jlist = firstneigh[i]; + const int jnum = numneigh[i]; + + int ninside = 0; + for (int jj = 0; jj < jnum; jj++) { + int j = jlist[jj]; + j &= NEIGHMASK; + const double delx = x[j][0] - xtmp; + const double dely = x[j][1] - ytmp; + const double delz = x[j][2] - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type[j]; + const int jelem = map[jtype]; + + if (rsq < descriptor->cutsq[ielem][jelem]) { + jatomdesc[ij] = j; + jelemdesc[ij] = jelem; + ij++; + } + } + } + + // compute descriptors + + descriptor->compute_descriptors(map, list, descriptors); + + // calculate descriptor gradients + + // descriptor->compute_descriptor_gradients(iatomdesc, ielemdesc, numneighdesc, + // jatomdesc, jelemdesc, graddesc, numneighdesc); + descriptor->compute_descriptor_gradients(map, list, graddesc, numneighdesc); + + // calculate force gradients w.r.t. parameters + + model->compute_force_gradients(descriptors, numlistdesc, iatomdesc, ielemdesc, + numneighdesc, jatomdesc, jelemdesc, graddesc, + yoffset, zoffset, gradforce, egradient); + + } // accumulate descriptor gradient contributions to global array @@ -343,7 +467,7 @@ void ComputeMLIAP::compute_array() own & ghost atoms ------------------------------------------------------------------------- */ -void ComputeMLIAP::dbdotr_compute() + void ComputeMLIAP::dbdotr_compute() { double **x = atom->x; int irow0 = 1+ndims_force*natoms; diff --git a/src/MLIAP/compute_mliap.h b/src/MLIAP/compute_mliap.h index 74b5fc1f93..780b0d0f76 100644 --- a/src/MLIAP/compute_mliap.h +++ b/src/MLIAP/compute_mliap.h @@ -31,6 +31,7 @@ class ComputeMLIAP : public Compute { void init(); void init_list(int, class NeighList *); void compute_array(); + void compute_array_alt(); double memory_usage(); private: @@ -42,17 +43,29 @@ class ComputeMLIAP : public Compute { double **gradforce; int *map; // map types to [0,nelements) int nelements; - + int gradgradflag; // 1 for graddesc, 0 for gamma double** descriptors; // descriptors for all atoms in list int ndescriptors; // number of descriptors - int gamma_max; // number of atoms allocated for beta, descriptors - int nparams; // number of model paramters per element + int nparams; // number of model parameters 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 double* egradient; // energy gradient w.r.t. parameters + // data structures for descriptor neighbor list + // this is for neighbors strictly inside deescriptor cutoff + + int numlistdesc_max; // number of atoms allocated in neighlist + int *numneighdesc; // numbers of neighbors for descriptors + int *iatomdesc; // list of descriptor atoms + int *ielemdesc; // list of descriptor elements + int nneighdesc_max; // number of ij neighbors allocated in graddesc + int *jatomdesc; // list of descriptor neighbor atoms + int *jelemdesc; // list of descriptor neighbor elements + int *neighdesc; // list of descriptor neighbors + double ***graddesc; // gradient of descriptors w.r.t. rij + class MLIAPModel* model; class MLIAPDescriptor* descriptor; diff --git a/src/MLIAP/mliap_descriptor.h b/src/MLIAP/mliap_descriptor.h index 072a5b4729..1b3813e7ac 100644 --- a/src/MLIAP/mliap_descriptor.h +++ b/src/MLIAP/mliap_descriptor.h @@ -26,8 +26,7 @@ public: virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int)=0; 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**, - double**, int, int)=0; + virtual void compute_descriptor_gradients(int*, NeighList*, double***, int*)=0; virtual void init()=0; virtual double memory_usage()=0; diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index c2879e942a..37ece96d7c 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -369,9 +369,7 @@ void MLIAPDescriptorSNAP::compute_gradients(int *map, NeighList* list, ---------------------------------------------------------------------- */ 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) + double ***graddesc, int *numneighdesc) { int i,j,jnum,ninside; double delx,dely,delz,evdwl,rsq; @@ -387,6 +385,7 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list numneigh = list->numneigh; firstneigh = list->firstneigh; + int ij = 0; for (int ii = 0; ii < list->inum; ii++) { i = list->ilist[ii]; @@ -431,6 +430,9 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list ninside++; } } + if (ninside != numneighdesc[i]) + error->one(FLERR,fmt::format("Screwed up neighbor count, ninside = {} jnum = {} i = {}", + ninside, numneighdesc[i], i)); if (chemflag) snaptr->compute_ui(ninside, ielem); @@ -458,13 +460,11 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(int *map, NeighList* list // Accumulate dB_k^i/dRi, dB_k^i/dRj for (int k = 0; k < ndescriptors; k++) { - graddesc[i][k] = snaptr->dblist[k][0]; - graddesc[i][k] = snaptr->dblist[k][1]; - graddesc[i][k] = snaptr->dblist[k][2]; - graddesc[j][k] = -snaptr->dblist[k][0]; - graddesc[j][k] = -snaptr->dblist[k][1]; - graddesc[j][k] = -snaptr->dblist[k][2]; + graddesc[ij][k][0] = snaptr->dblist[k][0]; + graddesc[ij][k][1] = snaptr->dblist[k][1]; + graddesc[ij][k][2] = snaptr->dblist[k][2]; } + ij++; } } diff --git a/src/MLIAP/mliap_descriptor_snap.h b/src/MLIAP/mliap_descriptor_snap.h index a929e01149..f33f513b88 100644 --- a/src/MLIAP/mliap_descriptor_snap.h +++ b/src/MLIAP/mliap_descriptor_snap.h @@ -26,8 +26,7 @@ public: virtual void compute_forces(class PairMLIAP*, class NeighList*, double**, int); 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**, - double**, int, int); + virtual void compute_descriptor_gradients(int*, NeighList*, double***, int*); virtual void init(); virtual double memory_usage(); diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h index 57f6f76214..8c13d896d2 100644 --- a/src/MLIAP/mliap_model.h +++ b/src/MLIAP/mliap_model.h @@ -26,6 +26,8 @@ public: virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int)=0; virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*)=0; virtual int get_gamma_nnz()=0; + virtual void compute_force_gradients(double**, int, int*, int*, int*, int*, + int*, double***, int, int, double**, double*)=0; virtual void init(); virtual double memory_usage(); int nelements; // # of unique elements diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp index b8d06155ba..07a9970b7e 100644 --- a/src/MLIAP/mliap_model_linear.cpp +++ b/src/MLIAP/mliap_model_linear.cpp @@ -137,6 +137,7 @@ void MLIAPModelLinear::param_gradient(int *map, NeighList* list, } } + /* ---------------------------------------------------------------------- count the number of non-zero entries in gamma matrix ---------------------------------------------------------------------- */ @@ -146,3 +147,46 @@ int MLIAPModelLinear::get_gamma_nnz() int inz = ndescriptors; return inz; } + +void MLIAPModelLinear::compute_force_gradients(double **descriptors, int numlistdesc, int *iatomdesc, int *ielemdesc, int *numneighdesc, int *jatomdesc, + int *jelemdesc, double ***graddesc, int yoffset, int zoffset, double **gradforce, double *egradient) { + + + // zero out energy gradients + + for (int l = 0; l < nelements*nparams; l++) + egradient[l] = 0.0; + + int ij = 0; + for (int ii = 0; ii < numlistdesc; ii++) { + const int i = iatomdesc[ii]; + const int ielem = ielemdesc[ii]; + const int elemoffset = nparams*ielem; + + for (int jj = 0; jj < numneighdesc[ii]; jj++) { + const int j = jatomdesc[ij]; + const int jelem = ielemdesc[ij]; + + int l = elemoffset+1; + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { + gradforce[i][l] += graddesc[ij][icoeff][0]; + gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]; + gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]; + gradforce[j][l] -= graddesc[ij][icoeff][0]; + gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]; + gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]; + l++; + } + ij++; + } + + // gradient of energy of atom I w.r.t. parameters + + int 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_linear.h b/src/MLIAP/mliap_model_linear.h index db4cb26514..eac9f574f2 100644 --- a/src/MLIAP/mliap_model_linear.h +++ b/src/MLIAP/mliap_model_linear.h @@ -26,7 +26,8 @@ public: virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int); virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*); virtual int get_gamma_nnz(); - + virtual void compute_force_gradients(double**, int, int*, int*, int*, int*, + int*, double***, int, int, double**, double*); protected: }; diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp index f5aae1d89d..ea7df974cb 100644 --- a/src/MLIAP/mliap_model_quadratic.cpp +++ b/src/MLIAP/mliap_model_quadratic.cpp @@ -217,3 +217,80 @@ int MLIAPModelQuadratic::get_gamma_nnz() return inz; } +void MLIAPModelQuadratic::compute_force_gradients(double **descriptors, int numlistdesc, int *iatomdesc, int *ielemdesc, int *numneighdesc, int *jatomdesc, + int *jelemdesc, double ***graddesc, int yoffset, int zoffset, double **gradforce, double *egradient) { + + + // zero out energy gradients + + for (int l = 0; l < nelements*nparams; l++) + egradient[l] = 0.0; + + int ij = 0; + for (int ii = 0; ii < numlistdesc; ii++) { + const int i = iatomdesc[ii]; + const int ielem = ielemdesc[ii]; + const int elemoffset = nparams*ielem; + + for (int jj = 0; jj < numneighdesc[ii]; jj++) { + const int j = jatomdesc[ij]; + const int jelem = ielemdesc[ij]; + const int ij0 = ij; + + int l = elemoffset+1; + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { + gradforce[i][l] += graddesc[ij][icoeff][0]; + gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]; + gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]; + gradforce[j][l] -= graddesc[ij][icoeff][0]; + gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]; + gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]; + l++; + } + + // quadratic contributions + + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { + double bveci = descriptors[ii][icoeff]; + gradforce[i][l] += graddesc[ij][icoeff][0]*bveci; + gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]*bveci; + gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]*bveci; + gradforce[j][l] -= graddesc[ij][icoeff][0]*bveci; + gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]*bveci; + gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]*bveci; + l++; + for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { + double bvecj = descriptors[ii][jcoeff]; + gradforce[i][l] += graddesc[ij][icoeff][0]*bvecj + graddesc[ij][jcoeff][0]*bveci; + gradforce[i][l+yoffset] += graddesc[ij][icoeff][1]*bvecj + graddesc[ij][jcoeff][1]*bveci; + gradforce[i][l+zoffset] += graddesc[ij][icoeff][2]*bvecj + graddesc[ij][jcoeff][2]*bveci; + gradforce[j][l] -= graddesc[ij][icoeff][0]*bvecj + graddesc[ij][jcoeff][0]*bveci; + gradforce[j][l+yoffset] -= graddesc[ij][icoeff][1]*bvecj + graddesc[ij][jcoeff][1]*bveci; + gradforce[j][l+zoffset] -= graddesc[ij][icoeff][2]*bvecj + graddesc[ij][jcoeff][2]*bveci; + l++; + } + } + ij++; + } + + // gradient of energy of atom I w.r.t. parameters + + int 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; + for (int jcoeff = icoeff+1; jcoeff < ndescriptors; jcoeff++) { + double bvecj = descriptors[ii][jcoeff]; + egradient[l++] += bveci*bvecj; + } + } + + } + +} diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h index b597fc7664..08cd1aed46 100644 --- a/src/MLIAP/mliap_model_quadratic.h +++ b/src/MLIAP/mliap_model_quadratic.h @@ -26,6 +26,8 @@ public: virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int); virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*); virtual int get_gamma_nnz(); + virtual void compute_force_gradients(double**, int, int*, int*, int*, int*, + int*, double***, int, int, double**, double*); protected: }; diff --git a/src/SNAP/compute_snap.cpp b/src/SNAP/compute_snap.cpp index 7eb216a63d..613fc1a862 100644 --- a/src/SNAP/compute_snap.cpp +++ b/src/SNAP/compute_snap.cpp @@ -212,16 +212,19 @@ void ComputeSnap::init() // find compute for reference energy - int ipe = modify->find_compute("thermo_pe"); + std::string id_pe = std::string("thermo_pe"); + int ipe = modify->find_compute(id_pe); if (ipe == -1) error->all(FLERR,"compute thermo_pe does not exist."); c_pe = modify->compute[ipe]; // add compute for reference virial tensor - modify->add_compute("snap_press all pressure NULL virial"); + std::string id_virial = std::string("snap_press"); + std::string pcmd = id_virial + " all pressure NULL virial"; + modify->add_compute(pcmd); - int ivirial = modify->find_compute("snap_press"); + int ivirial = modify->find_compute(id_virial); if (ivirial == -1) error->all(FLERR,"compute snap_press does not exist."); c_virial = modify->compute[ivirial];