ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Generated working example, not quite correct

Browse Source
This commit is contained in:
Aidan Thompson
2020-07-01 15:03:45 -06:00
parent 9c7005a91c
commit aee3d9f2cd
9 changed files with 217 additions and 46 deletions
Download Patch File Download Diff File Expand all files Collapse all files

95
examples/mliap/in.mliap.snap.compute Normal file
View File

@ -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}

66
src/MLIAP/compute_mliap.cpp
View File

@ -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

67
src/MLIAP/mliap_descriptor_snap.cpp
View File

@ -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;
}
/* ---------------------------------------------------------------------- */

10
src/MLIAP/mliap_model.cpp
View File

@ -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);

1
src/MLIAP/mliap_model.h
View File

@ -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;

13
src/MLIAP/mliap_model_linear.cpp
View File

@ -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;

1
src/MLIAP/mliap_model_linear.h
View File

@ -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*);

9
src/MLIAP/mliap_model_quadratic.cpp
View File

@ -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(){}
/* ----------------------------------------------------------------------

1
src/MLIAP/mliap_model_quadratic.h
View File

@ -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*);