diff --git a/examples/mliap/in.mliap.snap.compute b/examples/mliap/in.mliap.snap.compute new file mode 100644 index 0000000000..813e568387 --- /dev/null +++ b/examples/mliap/in.mliap.snap.compute @@ -0,0 +1,95 @@ +# Demonstrate bispectrum computes + +# initialize simulation + +variable nsteps index 0 +variable nrep equal 1 +variable a equal 2.0 +units metal + +# generate the box and atom positions using a BCC lattice + +variable nx equal ${nrep} +variable ny equal ${nrep} +variable nz equal ${nrep} + +boundary p p p + +atom_modify map hash +lattice bcc $a +region box block 0 ${nx} 0 ${ny} 0 ${nz} +create_box 2 box +create_atoms 2 box + +mass * 180.88 + +displace_atoms all random 0.1 0.1 0.1 123456 + +# set up reference potential + +variable zblcutinner equal 4 +variable zblcutouter equal 4.8 +variable zblz equal 73 +pair_style zbl ${zblcutinner} ${zblcutouter} +pair_coeff * * ${zblz} ${zblz} + +# choose SNA parameters + +variable twojmax equal 2 +variable rcutfac equal 1.0 +variable rfac0 equal 0.99363 +variable rmin0 equal 0 +variable radelem1 equal 2.3 +variable radelem2 equal 2.0 +variable wj1 equal 1.0 +variable wj2 equal 0.96 +variable quadratic equal 0 +variable bzero equal 0 +variable switch equal 0 +variable snap_options string & +"${rcutfac} ${rfac0} ${twojmax} ${radelem1} ${radelem2} ${wj1} ${wj2} rmin0 ${rmin0} quadraticflag ${quadratic} bzeroflag ${bzero} switchflag ${switch}" + +# set up per-atom computes + +compute b all sna/atom ${snap_options} +compute vb all snav/atom ${snap_options} +compute db all snad/atom ${snap_options} + +# perform sums over atoms + +group snapgroup1 type 1 +group snapgroup2 type 2 +compute bsum1 snapgroup1 reduce sum c_b[*] +compute bsum2 snapgroup2 reduce sum c_b[*] +# fix bsum1 all ave/time 1 1 1 c_bsum1 file bsum1.dat mode vector +# fix bsum2 all ave/time 1 1 1 c_bsum2 file bsum2.dat mode vector +compute vbsum all reduce sum c_vb[*] +# fix vbsum all ave/time 1 1 1 c_vbsum file vbsum.dat mode vector +variable db_2_30 equal c_db[2][30] + +# set up compute snap generating global array + +compute snap all mliap descriptor sna compute.mliap.descriptor model linear 2 6 +fix snap all ave/time 1 1 1 c_snap[*] file compute.snap.dat mode vector + +thermo 100 + +# test output: 1: total potential energy +# 2: xy component of stress tensor +# 3: Sum(B_{000}^i, all i of type 2) +# 4: xy component of Sum(Sum(r_j*dB_{222}^i/dR[j]), all i of type 2), all j) +# 5: z component of -Sum(d(B_{222}^i)/dR[2]), all i of type 2) +# +# followed by 5 counterparts from compute snap + +thermo_style custom & + pe pxy c_bsum2[1] c_vbsum[60] v_db_2_30 & + c_snap[1][11] c_snap[13][11] c_snap[1][6] c_snap[13][10] c_snap[7][10] +thermo_modify norm no + +# dump mydump_db all custom 1000 dump_db id c_db[*] +# dump_modify mydump_db sort id + +# Run MD + +run ${nsteps} diff --git a/src/MLIAP/compute_mliap.cpp b/src/MLIAP/compute_mliap.cpp index e86dce732e..69c00fc94c 100644 --- a/src/MLIAP/compute_mliap.cpp +++ b/src/MLIAP/compute_mliap.cpp @@ -35,14 +35,13 @@ enum{SCALAR,VECTOR,ARRAY}; ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), list(NULL), mliap(NULL), - mliap_peratom(NULL), mliapall(NULL) + mliap_peratom(NULL), mliapall(NULL), map(NULL), + descriptors(NULL), gamma_row_index(NULL), gamma_col_index(NULL), + gamma(NULL), egradient(NULL), model(NULL), descriptor(NULL) { - array_flag = 1; extarray = 0; - int ntypes = atom->ntypes; - if (narg < 4) error->all(FLERR,"Illegal compute mliap command"); @@ -53,19 +52,23 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : // process keywords - int iarg = 0; + int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"model") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command"); if (strcmp(arg[iarg+1],"linear") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command"); - model = new MLIAPModelLinear(lmp,arg[iarg+2]); - iarg += 3; + if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command"); + int ntmp1 = atoi(arg[iarg+2]); + int ntmp2 = atoi(arg[iarg+3]); + model = new MLIAPModelLinear(lmp,ntmp1,ntmp2); + iarg += 4; } else if (strcmp(arg[iarg+1],"quadratic") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal compute mliap command"); - model = new MLIAPModelQuadratic(lmp,arg[iarg+2]); - iarg += 3; + if (iarg+4 > narg) error->all(FLERR,"Illegal compute mliap command"); + int ntmp1 = atoi(arg[iarg+2]); + int ntmp2 = atoi(arg[iarg+3]); + model = new MLIAPModelQuadratic(lmp,ntmp1,ntmp2); + iarg += 4; } else error->all(FLERR,"Illegal compute mliap command"); modelflag = 1; } else if (strcmp(arg[iarg],"descriptor") == 0) { @@ -83,6 +86,7 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : if (modelflag == 0 || descriptorflag == 0) error->all(FLERR,"Illegal compute_style command"); + ndescriptors = descriptor->ndescriptors; nparams = model->nparams; nelements = model->nelements; gamma_nnz = model->get_gamma_nnz(); @@ -100,6 +104,8 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : size_peratom = ndims_peratom*nperdim*atom->ntypes; nmax = 0; + gamma_max = 0; + } /* ---------------------------------------------------------------------- */ @@ -147,14 +153,32 @@ void ComputeMLIAP::init() if (count > 1 && comm->me == 0) error->warning(FLERR,"More than one compute mliap"); + // initialize model and descriptor + + model->init(); + descriptor->init(); + + // consistency checks + + if (descriptor->ndescriptors != model->ndescriptors) + error->all(FLERR,"Incompatible model and descriptor definitions"); + if (descriptor->nelements != model->nelements) + error->all(FLERR,"Incompatible model and descriptor definitions"); + if (nelements != atom->ntypes) + error->all(FLERR,"nelements must equal ntypes"); + // allocate memory for global array + printf("size_array_rows = %d size_array_cols = %d\n",size_array_rows,size_array_cols); memory->create(mliap,size_array_rows,size_array_cols, "mliap:mliap"); memory->create(mliapall,size_array_rows,size_array_cols, "mliap:mliapall"); array = mliapall; + printf("nelements = %d nparams = %d\n",nelements,nparams); + memory->create(egradient,nelements*nparams,"ComputeMLIAP:egradient"); + // find compute for reference energy char *id_pe = (char *) "thermo_pe"; @@ -203,19 +227,11 @@ void ComputeMLIAP::compute_array() if (atom->nmax > nmax) { memory->destroy(mliap_peratom); nmax = atom->nmax; + printf("nmax = %d size_peratom = %d\n",nmax,size_peratom); memory->create(mliap_peratom,nmax,size_peratom, "mliap:mliap_peratom"); } - 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"); - memory->grow(egradient,nelements*nparams,"ComputeMLIAP:egradient"); - gamma_max = list->inum; - } - // clear global array for (int irow = 0; irow < size_array_rows; irow++) @@ -233,6 +249,15 @@ 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; + } + printf("gamma_max %d %d %d %d %d %d %p\n",gamma_max, ndescriptors, gamma_nnz, nelements, nparams, list->inum, list); + // compute descriptors, if needed if (model->nonlinearflag) @@ -312,6 +337,7 @@ void ComputeMLIAP::compute_array() int irow = 0; double reference_energy = c_pe->compute_scalar(); mliapall[irow++][lastcol] = reference_energy; + printf("Reference energy = %g %g %g %d %d\n",reference_energy,mliapall[irow-1][lastcol],array[irow-1][lastcol],irow-1,lastcol); // assign virial stress to last column // switch to Voigt notation diff --git a/src/MLIAP/mliap_descriptor_snap.cpp b/src/MLIAP/mliap_descriptor_snap.cpp index de3d7d3c9b..ffb4ab65ff 100644 --- a/src/MLIAP/mliap_descriptor_snap.cpp +++ b/src/MLIAP/mliap_descriptor_snap.cpp @@ -45,6 +45,13 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp, char *paramfilename): wjelem = NULL; snaptr = NULL; read_paramfile(paramfilename); + + snaptr = new SNA(lmp, rfac0, twojmax, + rmin0, switchflag, bzeroflag, + chemflag, bnormflag, wselfallflag, nelements); + + ndescriptors = snaptr->ncoeff; + } /* ---------------------------------------------------------------------- */ @@ -85,8 +92,9 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor const double ytmp = x[i][1]; const double ztmp = x[i][2]; const int itype = type[i]; - const int ielem = map[itype]; - const double radi = radelem[ielem]; + // element map not yet implemented + // const int ielem = map[itype]; + const int ielem = itype-1; jlist = list->firstneigh[i]; jnum = list->numneigh[i]; @@ -110,16 +118,18 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor delz = x[j][2] - ztmp; rsq = delx*delx + dely*dely + delz*delz; int jtype = type[j]; - int jelem = map[jtype]; + // element map not yet implemented + // const int jelem = map[jtype]; + const int jelem = jtype-1; if (rsq < cutsq[ielem][jelem]) { snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][1] = dely; snaptr->rij[ninside][2] = delz; snaptr->inside[ninside] = j; - snaptr->wj[ninside] = wjelem[jelem]; - snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); - snaptr->element[ninside] = jelem; // element index for chem snap + snaptr->wj[ninside] = wjelem[jelem]; + snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); + snaptr->element[ninside] = jelem; // element index for chem snap ninside++; } } @@ -129,6 +139,7 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor else snaptr->compute_ui(ninside, 0); snaptr->compute_zi(); + if (chemflag) snaptr->compute_bi(ielem); else @@ -168,7 +179,6 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double const double ztmp = x[i][2]; const int itype = type[i]; const int ielem = pairmliap->map[itype]; - const double radi = radelem[ielem]; jlist = firstneigh[i]; jnum = numneigh[i]; @@ -194,7 +204,7 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double int jtype = type[j]; int jelem = pairmliap->map[jtype]; - if (rsq < cutsq[itype][jtype]) { + if (rsq < cutsq[ielem][jelem]) { snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][1] = dely; snaptr->rij[ninside][2] = delz; @@ -286,8 +296,9 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list, const int itype = type[i]; int ielem = 0; if (chemflag) - ielem = map[itype]; - const double radi = radelem[ielem]; + // element map not yet implemented + // ielem = map[itype]; + const int ielem = itype-1; jlist = firstneigh[i]; jnum = numneigh[i]; @@ -311,28 +322,43 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list, delz = x[j][2] - ztmp; rsq = delx*delx + dely*dely + delz*delz; int jtype = type[j]; - int jelem = map[jtype]; - - if (rsq < cutsq[itype][jtype]) { + // element map not yet implemented + // int jelem = map[jtype]; + const int jelem = jtype-1; + if (rsq < cutsq[ielem][jelem]) { snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][1] = dely; snaptr->rij[ninside][2] = delz; snaptr->inside[ninside] = j; - snaptr->wj[ninside] = wjelem[jelem]; - snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); + snaptr->wj[ninside] = wjelem[jelem]; + snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); snaptr->element[ninside] = jelem; // element index for chem snap ninside++; } } - snaptr->compute_ui(ninside, ielem); + if(chemflag) + snaptr->compute_ui(ninside, ielem); + else + snaptr->compute_ui(ninside, 0); + + snaptr->compute_zi(); - snaptr->compute_bi(ielem); + if(chemflag) + snaptr->compute_bi(ielem); + else + snaptr->compute_bi(0); for (int jj = 0; jj < ninside; jj++) { const int j = snaptr->inside[jj]; + + if(chemflag) snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], snaptr->rcutij[jj], jj, snaptr->element[jj]); + else + snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], + snaptr->rcutij[jj], jj, 0); + snaptr->compute_dbidrj(); // Accumulate gamma_lk*dB_k/dRi, -gamma_lk**dB_k/dRj @@ -359,14 +385,7 @@ void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list, void MLIAPDescriptorSNAP::init() { - - snaptr = new SNA(lmp, rfac0, twojmax, - rmin0, switchflag, bzeroflag, - chemflag, bnormflag, wselfallflag, nelements); - snaptr->init(); - - ndescriptors = snaptr->ncoeff; } /* ---------------------------------------------------------------------- */ diff --git a/src/MLIAP/mliap_model.cpp b/src/MLIAP/mliap_model.cpp index 94aafc93a8..5ab96e19d4 100644 --- a/src/MLIAP/mliap_model.cpp +++ b/src/MLIAP/mliap_model.cpp @@ -40,6 +40,16 @@ MLIAPModel::MLIAPModel(LAMMPS* lmp, char* coefffilename) : Pointers(lmp) /* ---------------------------------------------------------------------- */ +MLIAPModel::MLIAPModel(LAMMPS* lmp, int nelements_in, int nparams_in) : Pointers(lmp) +{ + nelements = nelements_in; + nparams = nparams_in; + coeffelem = NULL; + nonlinearflag = 0; +} + +/* ---------------------------------------------------------------------- */ + MLIAPModel::~MLIAPModel() { memory->destroy(coeffelem); diff --git a/src/MLIAP/mliap_model.h b/src/MLIAP/mliap_model.h index 06b23707a9..57f6f76214 100644 --- a/src/MLIAP/mliap_model.h +++ b/src/MLIAP/mliap_model.h @@ -21,6 +21,7 @@ namespace LAMMPS_NS { class MLIAPModel : protected Pointers { public: MLIAPModel(LAMMPS*, char*); + MLIAPModel(LAMMPS*, int, int); ~MLIAPModel(); virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int)=0; virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*)=0; diff --git a/src/MLIAP/mliap_model_linear.cpp b/src/MLIAP/mliap_model_linear.cpp index 63e8220d14..a6d098ed26 100644 --- a/src/MLIAP/mliap_model_linear.cpp +++ b/src/MLIAP/mliap_model_linear.cpp @@ -31,7 +31,14 @@ using namespace LAMMPS_NS; MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, char* coefffilename) : MLIAPModel(lmp, coefffilename) { - nonlinearflag = 0; + ndescriptors = nparams - 1; +} + +/* ---------------------------------------------------------------------- */ + +MLIAPModelLinear::MLIAPModelLinear(LAMMPS* lmp, int nelements_in, int nparams_in) : + MLIAPModel(lmp, nelements_in, nparams_in) +{ ndescriptors = nparams - 1; } @@ -110,7 +117,9 @@ void MLIAPModelLinear::param_gradient(int *map, NeighList* list, i = list->ilist[ii]; const int itype = type[i]; - const int ielem = map[itype]; + // element map not yet implemented + // const int ielem = map[itype]; + const int ielem = itype-1; const int elemoffset = nparams*ielem; int l = elemoffset+1; diff --git a/src/MLIAP/mliap_model_linear.h b/src/MLIAP/mliap_model_linear.h index b7646aab45..db4cb26514 100644 --- a/src/MLIAP/mliap_model_linear.h +++ b/src/MLIAP/mliap_model_linear.h @@ -21,6 +21,7 @@ namespace LAMMPS_NS { class MLIAPModelLinear : public MLIAPModel { public: MLIAPModelLinear(LAMMPS*, char*); + MLIAPModelLinear(LAMMPS*, int, int); ~MLIAPModelLinear(); virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int); virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*); diff --git a/src/MLIAP/mliap_model_quadratic.cpp b/src/MLIAP/mliap_model_quadratic.cpp index 554d1a06a7..062e66959a 100644 --- a/src/MLIAP/mliap_model_quadratic.cpp +++ b/src/MLIAP/mliap_model_quadratic.cpp @@ -37,6 +37,15 @@ MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, char* coefffilename) : /* ---------------------------------------------------------------------- */ +MLIAPModelQuadratic::MLIAPModelQuadratic(LAMMPS* lmp, int nelements_in, int nparams_in) : + MLIAPModel(lmp, nelements_in, nparams_in) +{ + nonlinearflag = 1; + ndescriptors = sqrt(2*nparams)-1; +} + +/* ---------------------------------------------------------------------- */ + MLIAPModelQuadratic::~MLIAPModelQuadratic(){} /* ---------------------------------------------------------------------- diff --git a/src/MLIAP/mliap_model_quadratic.h b/src/MLIAP/mliap_model_quadratic.h index 4cfc669bbb..b597fc7664 100644 --- a/src/MLIAP/mliap_model_quadratic.h +++ b/src/MLIAP/mliap_model_quadratic.h @@ -21,6 +21,7 @@ namespace LAMMPS_NS { class MLIAPModelQuadratic : public MLIAPModel { public: MLIAPModelQuadratic(LAMMPS*, char*); + MLIAPModelQuadratic(LAMMPS*, int, int); ~MLIAPModelQuadratic(); virtual void gradient(class PairMLIAP*, class NeighList*, double**, double**, int); virtual void param_gradient(int*, class NeighList*, double**, int**, int**, double**, double*);