Commit 9c7005a9 authored by Aidan Thompson's avatar Aidan Thompson
Browse files

It compiles, but not yet working

parent cfa12bc1

src/MLIAP/compute_mliap.cpp

+42 −18
Original line number Diff line number Diff line
@@ -34,7 +34,7 @@ using namespace LAMMPS_NS; 
enum{SCALAR,VECTOR,ARRAY};



ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) :

  Compute(lmp, narg, arg), cutsq(NULL), list(NULL), mliap(NULL),

  Compute(lmp, narg, arg), list(NULL), mliap(NULL),

  mliap_peratom(NULL), mliapall(NULL)

{
@@ -46,6 +46,11 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : 
  if (narg < 4)

    error->all(FLERR,"Illegal compute mliap command");



  // set flags for required keywords



  int modelflag = 0;

  int descriptorflag = 0;



  // process keywords



  int iarg = 0;
@@ -62,6 +67,7 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : 
        model = new MLIAPModelQuadratic(lmp,arg[iarg+2]);

        iarg += 3;

      } else error->all(FLERR,"Illegal compute mliap command");

      modelflag = 1;

    } else if (strcmp(arg[iarg],"descriptor") == 0) {

      if (iarg+2 > narg) error->all(FLERR,"Illegal compute mliap command");

      if (strcmp(arg[iarg+1],"sna") == 0) {
@@ -69,10 +75,17 @@ ComputeMLIAP::ComputeMLIAP(LAMMPS *lmp, int narg, char **arg) : 
        descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]);

        iarg += 3;

      } else error->all(FLERR,"Illegal compute mliap command");

      descriptorflag = 1;

    } else

      error->all(FLERR,"Illegal compute mliap command");

  }

  }



  if (modelflag == 0 || descriptorflag == 0)

    error->all(FLERR,"Illegal compute_style command");



  nparams = model->nparams;

  nelements = model->nelements;

  gamma_nnz = model->get_gamma_nnz();

  nperdim = nparams;

  ndims_force = 3;

  ndims_virial = 6;
@@ -96,9 +109,17 @@ ComputeMLIAP::~ComputeMLIAP() 
  memory->destroy(mliap);

  memory->destroy(mliapall);

  memory->destroy(mliap_peratom);

  memory->destroy(cutsq);



  memory->destroy(map);



  memory->destroy(descriptors);

  memory->destroy(gamma_row_index);

  memory->destroy(gamma_col_index);

  memory->destroy(gamma);

  memory->destroy(egradient);



  delete model;

  delete descriptor;

}



/* ---------------------------------------------------------------------- */
@@ -108,7 +129,7 @@ void ComputeMLIAP::init() 
  if (force->pair == NULL)

    error->all(FLERR,"Compute mliap requires a pair style be defined");



  if (descriptor->get_cutmax() > force->pair->cutforce)

  if (descriptor->cutmax > force->pair->cutforce)

    error->all(FLERR,"Compute mliap cutoff is longer than pairwise cutoff");



  // need an occasional full neighbor list
@@ -187,8 +208,11 @@ void ComputeMLIAP::compute_array() 
  }



  if (gamma_max < list->inum) {

    memory->grow(descriptors,list->inum,ndescriptors,"PairMLIAP:descriptors");

    memory->grow(gamma,nparams,list->inum,ndescriptors,"PairMLIAP:gamma");

    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;

  }
@@ -214,19 +238,12 @@ void ComputeMLIAP::compute_array() 
  if (model->nonlinearflag)

    descriptor->forward(map, list, descriptors);



  // ***********THIS IS NOT RIGHT**********************

  // This whole idea is flawed. The gamma matrix is too big to

  // store. Instead, we should generate the A matrix,

  // just as ComputeSNAP does, and then pass it to 

  // the model, which can evaluate gradients of E, F, sigma,

  // w.r.t. model parameters.



  // calculate descriptor contributions to parameter gradients

  // and gamma = double gradient w.r.t. parameters and descriptors



  // i.e. gamma = d2E/d\sigma.dB_i

  // sigma is a parameter and B_i is a descriptor of atom i

  // for SNAP, this is a sparse nparams*natoms*ndescriptors matrix,

  // 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. 

 

  // *******Not implemented yet*****************
@@ -237,16 +254,23 @@ void ComputeMLIAP::compute_array() 
  // For the quadratic model, the energy row will be similar,

  // while gamma will be 1's, 0's and Bi's  



  //   model->param_gradient(list, descriptors, mliap[0], gamma);

  model->param_gradient(map, list, descriptors, gamma_row_index, 

                        gamma_col_index, gamma, egradient);





  // calculate descriptor gradient contributions to parameter gradients



  // *******Not implemented yet*****************

  // This will just take gamma and multiply it with

  // descriptor gradient contributions i.e. dblist

  // this will resemble snadi accumualation in ComputeSNAP

  // this will resemble snadi accumulation in ComputeSNAP



  // descriptor->param_backward(list, gamma, snadi);



  descriptor->param_backward(map, list, gamma_nnz, gamma_row_index, 

                             gamma_col_index, gamma, mliap_peratom,

                             yoffset, zoffset);



  // accumulate descriptor gradient contributions to global array



  for (int itype = 0; itype < atom->ntypes; itype++) {

src/MLIAP/compute_mliap.h

+5 −2
Original line number Diff line number Diff line
@@ -37,18 +37,21 @@ class ComputeMLIAP : public Compute { 
  int natoms, nmax, size_peratom, lastcol;

  int nperdim, yoffset, zoffset;

  int ndims_peratom, ndims_force, ndims_virial;

  double **cutsq;

  class NeighList *list;

  double **mliap, **mliapall;

  double **mliap_peratom;

  int *map;  // map types to [0,nelements)

  int nelements;



  double*** gamma;             // gammas for all atoms in list

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



  class MLIAPModel* model;

  class MLIAPDescriptor* descriptor;

src/MLIAP/mliap_descriptor.h

+4 −3
Original line number Diff line number Diff line
@@ -24,15 +24,16 @@ public: 
  ~MLIAPDescriptor();

  virtual void forward(int*, class NeighList*, double**)=0;

  virtual void backward(class PairMLIAP*, class NeighList*, double**, int)=0;

  virtual void param_backward(int*, class NeighList*, int, int**, int**, double**, 

                              double**, int, int)=0;

  virtual void init()=0;

  virtual double get_cutoff(int, int)=0;

  virtual double get_cutmax()=0;

  virtual double memory_usage()=0;



  int ndescriptors;              // number of descriptors

  int nelements;                 // # of unique elements

  char **elements;               // names of unique elements



  double **cutsq;                // nelem x nelem rcutsq values 

  double cutmax;                 // maximum cutoff needed

protected:



};

src/MLIAP/mliap_descriptor_snap.cpp

+116 −29
Original line number Diff line number Diff line
@@ -58,6 +58,7 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() 
    delete[] elements;

    memory->destroy(radelem);

    memory->destroy(wjelem);

    memory->destroy(cutsq);

  }



  delete snaptr;
@@ -111,16 +112,13 @@ void MLIAPDescriptorSNAP::forward(int* map, NeighList* list, double **descriptor 
      int jtype = type[j];

      int jelem = map[jtype];



      //      printf("i = %d j = %d itype = %d jtype = %d cutsq[i][j] = %g rsq = %g\n",i,j,itype,jtype,cutsq[itype][jtype],rsq);



      double rcutsqtmp = get_cutoff(ielem, jelem);

      if (rsq < rcutsqtmp*rcutsqtmp) {

      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] = (radi + radelem[jelem])*rcutfac;

        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);

        snaptr->element[ninside] = jelem; // element index for chem snap

        ninside++;

      }
@@ -196,13 +194,13 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double 
      int jtype = type[j];

      int jelem = pairmliap->map[jtype];



      if (rsq < pairmliap->cutsq[itype][jtype]&&rsq>1e-20) {

      if (rsq < cutsq[itype][jtype]) {

        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] = (radi + radelem[jelem])*rcutfac;

        snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]);

        snaptr->element[ninside] = jelem; // element index for chem snap

        ninside++;

      }
@@ -242,7 +240,7 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double 
      f[j][1] -= fij[1];

      f[j][2] -= fij[2];



      // add in gloabl and per-atom virial contributions

      // add in global and per-atom virial contributions

      // this is optional and has no effect on force calculation

      

      if (vflag)
@@ -256,6 +254,105 @@ void MLIAPDescriptorSNAP::backward(PairMLIAP* pairmliap, NeighList* list, double 


}



/* ----------------------------------------------------------------------

   compute forces for each atom

   ---------------------------------------------------------------------- */



void MLIAPDescriptorSNAP::param_backward(int *map, NeighList* list, 

                                         int gamma_nnz, int **gamma_row_index, 

                                         int **gamma_col_index, double **gamma, double **snadi,

                                         int yoffset, int zoffset)

{

  int i,j,jnum,ninside;

  double delx,dely,delz,evdwl,rsq;

  double fij[3];

  int *jlist,*numneigh,**firstneigh;



  double **x = atom->x;

  double **f = atom->f;

  int *type = atom->type;

  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  numneigh = list->numneigh;

  firstneigh = list->firstneigh;



  for (int ii = 0; ii < list->inum; ii++) {

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

    int ielem = 0;

    if (chemflag)

      ielem = map[itype];

    const double radi = radelem[ielem];



    jlist = firstneigh[i];

    jnum = numneigh[i];



    // insure rij, inside, wj, and rcutij are of size jnum



    snaptr->grow_rij(jnum);



    // rij[][3] = displacements between atom I and those neighbors

    // inside = indices of neighbors of I within cutoff

    // wj = weights for neighbors of I within cutoff

    // rcutij = cutoffs for neighbors of I within cutoff

    // note Rij sign convention => dU/dRij = dU/dRj = -dU/dRi



    ninside = 0;

    for (int jj = 0; jj < jnum; jj++) {

      j = jlist[jj];

      j &= NEIGHMASK;

      delx = x[j][0] - xtmp;

      dely = x[j][1] - ytmp;

      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]) {

        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

        ninside++;

      }

    }



    snaptr->compute_ui(ninside, ielem);

    snaptr->compute_zi();

    snaptr->compute_bi(ielem);



    for (int jj = 0; jj < ninside; jj++) {

      const int j = snaptr->inside[jj];

      snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],

                             snaptr->rcutij[jj], jj, snaptr->element[jj]);

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

          snadi[i][l]         += gamma[ii][inz]*snaptr->dblist[k][0];

          snadi[i][l+yoffset] += gamma[ii][inz]*snaptr->dblist[k][1];

          snadi[i][l+zoffset] += gamma[ii][inz]*snaptr->dblist[k][2];

          snadi[j][l]         -= gamma[ii][inz]*snaptr->dblist[k][0];

          snadi[j][l+yoffset] -= gamma[ii][inz]*snaptr->dblist[k][1];

          snadi[j][l+zoffset] -= gamma[ii][inz]*snaptr->dblist[k][2];

        }



    }

  }



}



/* ----------------------------------------------------------------------

   set coeffs for one or more type pairs

------------------------------------------------------------------------- */
@@ -270,7 +367,6 @@ void MLIAPDescriptorSNAP::init() 
  snaptr->init();



  ndescriptors = snaptr->ncoeff;



}



/* ---------------------------------------------------------------------- */
@@ -417,31 +513,22 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) 
      !elementsflag || !radelemflag || !wjelemflag)

    error->all(FLERR,"Incorrect SNAP parameter file");



}

  // construct cutsq



/* ----------------------------------------------------------------------

   provide cutoff distance for two elements

------------------------------------------------------------------------- */



double MLIAPDescriptorSNAP::get_cutoff(int ielem, int jelem)

{

  return (radelem[ielem] + radelem[jelem])*rcutfac;

}



/* ----------------------------------------------------------------------

   calculate maximum cutoff distance

------------------------------------------------------------------------- */



double MLIAPDescriptorSNAP::get_cutmax()

{

  double cut;

  double cutmax = 0.0;

  for(int ielem = 0; ielem <= nelements; ielem++) {

  cutmax = 0.0;

  memory->create(cutsq,nelements,nelements,"mliap/descriptor/snap:cutsq");

  for (int ielem = 0; ielem < nelements; ielem++) {

    cut = 2.0*radelem[ielem]*rcutfac;

    if (cut > cutmax) cutmax = cut;

    return cutmax;

    cutsq[ielem][ielem] = cut*cut;

    printf("ielem = %d ielem = %d cutsq[i][i] = %g\n",ielem, ielem, cutsq[ielem][ielem]);

    for(int jelem = ielem+1; jelem < nelements; jelem++) {

      cut = (radelem[ielem]+radelem[jelem])*rcutfac;

      cutsq[ielem][jelem] = cutsq[jelem][ielem] = cut*cut;

      printf("ielem = %d jelem = %d cutsq[i][j] = %g\n",ielem, jelem, cutsq[ielem][jelem]);

    }

  }

  return cutmax;

}



/* ----------------------------------------------------------------------

src/MLIAP/mliap_descriptor_snap.h

+2 −2
Original line number Diff line number Diff line
@@ -24,9 +24,9 @@ public: 
  ~MLIAPDescriptorSNAP();

  virtual void forward(int*, class NeighList*, double**);

  virtual void backward(class PairMLIAP*, class NeighList*, double**, int);

  virtual void param_backward(int*, class NeighList*, int, int**, int**, double**, 

                              double**, int, int);

  virtual void init();

  virtual double get_cutoff(int, int);

  virtual double get_cutmax();

  virtual double memory_usage();



  double rcutfac;                // declared public to workaround gcc 4.9

CRA Git | Maintained and supported by SUSTech CRA and CCSE