Commit 06a199f7 authored by Aidan Thompson's avatar Aidan Thompson
Browse files

Added gradgradflag option to compute mliap

parent 6bf32909

examples/mliap/in.mliap.quadratic.compute

+1 −1
Original line number Diff line number Diff line
@@ -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

examples/mliap/in.mliap.snap.compute

+1 −1
Original line number Diff line number Diff line
@@ -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

src/MLIAP/compute_mliap.cpp

+154 −30
Original line number Diff line number Diff line
@@ -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,14 +265,20 @@ 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;

  }



  if (gradgradflag) {



    // compute descriptors

    

    descriptor->compute_descriptors(map, list, descriptors);
@@ -275,6 +301,104 @@ void ComputeMLIAP::compute_array() 
                                  gamma_col_index, gamma, gradforce,

                                  yoffset, zoffset);

    

  } else {



    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



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

src/MLIAP/compute_mliap.h

+16 −3
Original line number Diff line number Diff line
@@ -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;

src/MLIAP/mliap_descriptor.h

+1 −2
Original line number Diff line number Diff line
@@ -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;

CRA Git | Maintained and supported by SUSTech CRA and CCSE