Commit 3dd9ad48 authored by Aidan Thompson's avatar Aidan Thompson
Browse files

Turned on MLIAPModelSNAP

parent 8789ff65

src/MLIAP/mliap_descriptor_snap.cpp

0 → 100644
+751 −0
Original line number Diff line number Diff line 
/* ----------------------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov



   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.



   See the README file in the top-level LAMMPS directory.

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



#include "mliap_descriptor_snap.h"

#include <mpi.h>

#include <cmath>

#include <cstdlib>

#include <cstring>

#include "atom.h"

#include "force.h"

#include "comm.h"

#include "neighbor.h"

#include "neigh_list.h"

#include "neigh_request.h"

#include "sna.h"

#include "memory.h"

#include "error.h"



using namespace LAMMPS_NS;



#define MAXLINE 1024

#define MAXWORD 3



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



//MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp): Pointers(lmp)

MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *lmp): Pair(lmp)

{

  single_enable = 0;

  restartinfo = 0;

  one_coeff = 1;

  manybody_flag = 1;



  nelements = 0;

  elements = NULL;

  radelem = NULL;

  wjelem = NULL;

  coeffelem = NULL;



  beta_max = 0;

  beta = NULL;

  bispectrum = NULL;

  snaptr = NULL;

}



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



MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP()

{

  if (copymode) return;



  if (nelements) {

    for (int i = 0; i < nelements; i++)

      delete[] elements[i];

    delete[] elements;

    memory->destroy(radelem);

    memory->destroy(wjelem);

    memory->destroy(coeffelem);

  }



  memory->destroy(beta);

  memory->destroy(bispectrum);



  delete snaptr;



  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(cutsq);

    memory->destroy(map);

  }



}



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

   This version is a straightforward implementation

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



void MLIAPDescriptorSNAP::compute(int eflag, int vflag)

{

  int i,j,jnum,ninside;

  double delx,dely,delz,evdwl,rsq;

  double fij[3];

  int *jlist,*numneigh,**firstneigh;



  ev_init(eflag,vflag);



  double **x = atom->x;

  double **f = atom->f;

  int *type = atom->type;

  int nlocal = atom->nlocal;

  int newton_pair = force->newton_pair;



  if (beta_max < list->inum) {

    memory->grow(beta,list->inum,ncoeff,"MLIAPDescriptorSNAP:beta");

    memory->grow(bispectrum,list->inum,ncoeff,"MLIAPDescriptorSNAP:bispectrum");

    beta_max = list->inum;

  }



  // compute dE_i/dB_i = beta_i for all i in list



  if (quadraticflag || eflag)

    compute_bispectrum();

  compute_beta();



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

    const int 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]&&rsq>1e-20) {

        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->element[ninside] = jelem; // element index for alloy snap

        ninside++;

      }

    }



    // compute Ui, Yi for atom I



    snaptr->compute_ui(ninside, ielem);



    // for neighbors of I within cutoff:

    // compute Fij = dEi/dRj = -dEi/dRi

    // add to Fi, subtract from Fj



    snaptr->compute_yi(beta[ii]);

    //for (int q=0; q<snaptr->idxu_max*2; q++){

    //  fprintf(screen, "%i %f\n",q, snaptr->ylist_r[q]);

    //}



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

      int j = snaptr->inside[jj];

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

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



      snaptr->compute_deidrj(fij);



      f[i][0] += fij[0];

      f[i][1] += fij[1];

      f[i][2] += fij[2];

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

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

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



      // tally per-atom virial contribution



      if (vflag)

        ev_tally_xyz(i,j,nlocal,newton_pair,0.0,0.0,

                     fij[0],fij[1],fij[2],

                     -snaptr->rij[jj][0],-snaptr->rij[jj][1],

                     -snaptr->rij[jj][2]);

    }

    //fprintf(screen, "%f %f %f\n",f[i][0],f[i][1],f[i][2]);



    // tally energy contribution



    if (eflag) {



      // evdwl = energy of atom I, sum over coeffs_k * Bi_k



      double* coeffi = coeffelem[ielem];

      evdwl = coeffi[0];

      // snaptr->copy_bi2bvec();



      // E = beta.B + 0.5*B^t.alpha.B



      // linear contributions



      for (int icoeff = 0; icoeff < ncoeff; icoeff++)

        evdwl += coeffi[icoeff+1]*bispectrum[ii][icoeff];



      // quadratic contributions



      if (quadraticflag) {

        int k = ncoeff+1;

        for (int icoeff = 0; icoeff < ncoeff; icoeff++) {

          double bveci = bispectrum[ii][icoeff];

          evdwl += 0.5*coeffi[k++]*bveci*bveci;

          for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {

            double bvecj = bispectrum[ii][jcoeff];

            evdwl += coeffi[k++]*bveci*bvecj;

          }

        }

      }

      ev_tally_full(i,2.0*evdwl,0.0,0.0,0.0,0.0,0.0);

    }



  }



  if (vflag_fdotr) virial_fdotr_compute();

}



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

   compute beta

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



void MLIAPDescriptorSNAP::compute_beta()

{

  int i;

  int *type = atom->type;



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

    i = list->ilist[ii];

    const int itype = type[i];

    const int ielem = map[itype];

    double* coeffi = coeffelem[ielem];



    for (int icoeff = 0; icoeff < ncoeff; icoeff++)

      beta[ii][icoeff] = coeffi[icoeff+1];



    if (quadraticflag) {

      int k = ncoeff+1;

      for (int icoeff = 0; icoeff < ncoeff; icoeff++) {

        double bveci = bispectrum[ii][icoeff];

        beta[ii][icoeff] += coeffi[k]*bveci;

        k++;

        for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {

          double bvecj = bispectrum[ii][jcoeff];

          beta[ii][icoeff] += coeffi[k]*bvecj;

          beta[ii][jcoeff] += coeffi[k]*bveci;

          k++;

        }

      }

    }

  }

}



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

   compute bispectrum

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



void MLIAPDescriptorSNAP::compute_bispectrum()

{

  int i,j,jnum,ninside;

  double delx,dely,delz,rsq;

  int *jlist;



  double **x = atom->x;

  int *type = atom->type;



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

    const int ielem = map[itype];

    const double radi = radelem[ielem];



    jlist = list->firstneigh[i];

    jnum = list->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]&&rsq>1e-20) {

        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->element[ninside] = jelem; // element index for alloy snap

        ninside++;

      }

    }



    snaptr->compute_ui(ninside, ielem);

    snaptr->compute_zi();

    snaptr->compute_bi(ielem);



    for (int icoeff = 0; icoeff < ncoeff; icoeff++){

      bispectrum[ii][icoeff] = snaptr->blist[icoeff];

    }

  }



}



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

   allocate all arrays

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



void MLIAPDescriptorSNAP::allocate()

{

  allocated = 1;

  int n = atom->ntypes;



  memory->create(setflag,n+1,n+1,"mliap_descriptor_snap:setflag");

  memory->create(cutsq,n+1,n+1,"mliap_descriptor_snap:cutsq");

  memory->create(map,n+1,"mliap_descriptor_snap:map");

}



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

   global settings

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



void MLIAPDescriptorSNAP::settings(int narg, char ** /* arg */)

{

  if (narg > 0)

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

}



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

   set coeffs for one or more type pairs

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



void MLIAPDescriptorSNAP::coeff(int narg, char **arg)

{

  if (narg < 5) error->all(FLERR,"Incorrect args for pair coefficients");

  if (!allocated) allocate();



  if (nelements) {

    for (int i = 0; i < nelements; i++)

      delete[] elements[i];

    delete[] elements;

    memory->destroy(radelem);

    memory->destroy(wjelem);

    memory->destroy(coeffelem);

  }



  char* type1 = arg[0];

  char* type2 = arg[1];

  char* coefffilename = arg[2];

  char* paramfilename = arg[3];

  char** elemtypes = &arg[4];



  // insure I,J args are * *



  if (strcmp(type1,"*") != 0 || strcmp(type2,"*") != 0)

    error->all(FLERR,"Incorrect args for pair coefficients");



  // read snapcoeff and snapparam files



  read_files(coefffilename,paramfilename);



  if (!quadraticflag)

    ncoeff = ncoeffall - 1;

  else {



    // ncoeffall should be (ncoeff+2)*(ncoeff+1)/2

    // so, ncoeff = floor(sqrt(2*ncoeffall))-1



    ncoeff = sqrt(2*ncoeffall)-1;

    ncoeffq = (ncoeff*(ncoeff+1))/2;

    int ntmp = 1+ncoeff+ncoeffq;

    if (ntmp != ncoeffall) {

      printf("ncoeffall = %d ntmp = %d ncoeff = %d \n",ncoeffall,ntmp,ncoeff);

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

    }

  }



  // read args that map atom types to SNAP elements

  // map[i] = which element the Ith atom type is, -1 if not mapped

  // map[0] is not used



  for (int i = 1; i <= atom->ntypes; i++) {

    char* elemname = elemtypes[i-1];

    int jelem;

    for (jelem = 0; jelem < nelements; jelem++)

      if (strcmp(elemname,elements[jelem]) == 0)

        break;



    if (jelem < nelements)

      map[i] = jelem;

    else if (strcmp(elemname,"NULL") == 0) map[i] = -1;

    else error->all(FLERR,"Incorrect args for pair coefficients");

  }



  // clear setflag since coeff() called once with I,J = * *



  int n = atom->ntypes;

  for (int i = 1; i <= n; i++)

    for (int j = i; j <= n; j++)

      setflag[i][j] = 0;



  // set setflag i,j for type pairs where both are mapped to elements



  int count = 0;

  for (int i = 1; i <= n; i++)

    for (int j = i; j <= n; j++)

      if (map[i] >= 0 && map[j] >= 0) {

        setflag[i][j] = 1;

        count++;

      }



  if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");



  snaptr = new SNA(lmp, rfac0, twojmax,

                   rmin0, switchflag, bzeroflag,

                   alloyflag, wselfallflag, nelements);



  if (ncoeff != snaptr->ncoeff) {

    if (comm->me == 0)

      printf("ncoeff = %d snancoeff = %d \n",ncoeff,snaptr->ncoeff);

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

  }



  // Calculate maximum cutoff for all elements



  rcutmax = 0.0;

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

    rcutmax = MAX(2.0*radelem[ielem]*rcutfac,rcutmax);



}



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

   init specific to this pair style

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



void MLIAPDescriptorSNAP::init_style()

{

  if (force->newton_pair == 0)

    error->all(FLERR,"Pair style SNAP requires newton pair on");



  // need a full neighbor list



  int irequest = neighbor->request(this,instance_me);

  neighbor->requests[irequest]->half = 0;

  neighbor->requests[irequest]->full = 1;



  snaptr->init();



}



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

   init for one type pair i,j and corresponding j,i

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



double MLIAPDescriptorSNAP::init_one(int i, int j)

{

  if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");

  return (radelem[map[i]] +

          radelem[map[j]])*rcutfac;

}



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



void MLIAPDescriptorSNAP::read_files(char *coefffilename, char *paramfilename)

{



  // open SNAP coefficient file on proc 0



  FILE *fpcoeff;

  if (comm->me == 0) {

    fpcoeff = force->open_potential(coefffilename);

    if (fpcoeff == NULL) {

      char str[128];

      snprintf(str,128,"Cannot open SNAP coefficient file %s",coefffilename);

      error->one(FLERR,str);

    }

  }



  char line[MAXLINE],*ptr;

  int eof = 0;



  int n;

  int nwords = 0;

  while (nwords == 0) {

    if (comm->me == 0) {

      ptr = fgets(line,MAXLINE,fpcoeff);

      if (ptr == NULL) {

        eof = 1;

        fclose(fpcoeff);

      } else n = strlen(line) + 1;

    }

    MPI_Bcast(&eof,1,MPI_INT,0,world);

    if (eof) break;

    MPI_Bcast(&n,1,MPI_INT,0,world);

    MPI_Bcast(line,n,MPI_CHAR,0,world);



    // strip comment, skip line if blank



    if ((ptr = strchr(line,'#'))) *ptr = '\0';

    nwords = atom->count_words(line);

  }

  if (nwords != 2)

    error->all(FLERR,"Incorrect format in SNAP coefficient file");



  // words = ptrs to all words in line

  // strip single and double quotes from words



  char* words[MAXWORD];

  int iword = 0;

  words[iword] = strtok(line,"' \t\n\r\f");

  iword = 1;

  words[iword] = strtok(NULL,"' \t\n\r\f");



  nelements = atoi(words[0]);

  ncoeffall = atoi(words[1]);



  // set up element lists



  elements = new char*[nelements];

  memory->create(radelem,nelements,"mliap_snap_descriptor:radelem");

  memory->create(wjelem,nelements,"mliap_snap_descriptor:wjelem");

  memory->create(coeffelem,nelements,ncoeffall,"mliap_snap_descriptor:coeffelem");



  // Loop over nelements blocks in the SNAP coefficient file



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



    if (comm->me == 0) {

      ptr = fgets(line,MAXLINE,fpcoeff);

      if (ptr == NULL) {

        eof = 1;

        fclose(fpcoeff);

      } else n = strlen(line) + 1;

    }

    MPI_Bcast(&eof,1,MPI_INT,0,world);

    if (eof)

      error->all(FLERR,"Incorrect format in SNAP coefficient file");

    MPI_Bcast(&n,1,MPI_INT,0,world);

    MPI_Bcast(line,n,MPI_CHAR,0,world);



    nwords = atom->count_words(line);

    if (nwords != 3)

      error->all(FLERR,"Incorrect format in SNAP coefficient file");



    iword = 0;

    words[iword] = strtok(line,"' \t\n\r\f");

    iword = 1;

    words[iword] = strtok(NULL,"' \t\n\r\f");

    iword = 2;

    words[iword] = strtok(NULL,"' \t\n\r\f");



    char* elemtmp = words[0];

    int n = strlen(elemtmp) + 1;

    elements[ielem] = new char[n];

    strcpy(elements[ielem],elemtmp);



    radelem[ielem] = atof(words[1]);

    wjelem[ielem] = atof(words[2]);





    if (comm->me == 0) {

      if (screen) fprintf(screen,"SNAP Element = %s, Radius %g, Weight %g \n",

                          elements[ielem], radelem[ielem], wjelem[ielem]);

      if (logfile) fprintf(logfile,"SNAP Element = %s, Radius %g, Weight %g \n",

                          elements[ielem], radelem[ielem], wjelem[ielem]);

    }



    for (int icoeff = 0; icoeff < ncoeffall; icoeff++) {

      if (comm->me == 0) {

        ptr = fgets(line,MAXLINE,fpcoeff);

        if (ptr == NULL) {

          eof = 1;

          fclose(fpcoeff);

        } else n = strlen(line) + 1;

      }



      MPI_Bcast(&eof,1,MPI_INT,0,world);

      if (eof)

        error->all(FLERR,"Incorrect format in SNAP coefficient file");

      MPI_Bcast(&n,1,MPI_INT,0,world);

      MPI_Bcast(line,n,MPI_CHAR,0,world);



      nwords = atom->count_words(line);

      if (nwords != 1)

        error->all(FLERR,"Incorrect format in SNAP coefficient file");



      iword = 0;

      words[iword] = strtok(line,"' \t\n\r\f");



      coeffelem[ielem][icoeff] = atof(words[0]);



    }

  }



  if (comm->me == 0) fclose(fpcoeff);



  // set flags for required keywords



  rcutfacflag = 0;

  twojmaxflag = 0;



  // Set defaults for optional keywords



  rfac0 = 0.99363;

  rmin0 = 0.0;

  switchflag = 1;

  bzeroflag = 1;

  bnormflag = 0;

  quadraticflag = 0;

  alloyflag = 0;

  wselfallflag = 0;

  chunksize = 2000;



  // open SNAP parameter file on proc 0



  FILE *fpparam;

  if (comm->me == 0) {

    fpparam = force->open_potential(paramfilename);

    if (fpparam == NULL) {

      char str[128];

      snprintf(str,128,"Cannot open SNAP parameter file %s",paramfilename);

      error->one(FLERR,str);

    }

  }



  eof = 0;

  while (1) {

    if (comm->me == 0) {

      ptr = fgets(line,MAXLINE,fpparam);

      if (ptr == NULL) {

        eof = 1;

        fclose(fpparam);

      } else n = strlen(line) + 1;

    }

    MPI_Bcast(&eof,1,MPI_INT,0,world);

    if (eof) break;

    MPI_Bcast(&n,1,MPI_INT,0,world);

    MPI_Bcast(line,n,MPI_CHAR,0,world);



    // strip comment, skip line if blank



    if ((ptr = strchr(line,'#'))) *ptr = '\0';

    nwords = atom->count_words(line);

    if (nwords == 0) continue;



    if (nwords != 2)

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



    // words = ptrs to all words in line

    // strip single and double quotes from words



    char* keywd = strtok(line,"' \t\n\r\f");

    char* keyval = strtok(NULL,"' \t\n\r\f");



    if (comm->me == 0) {

      if (screen) fprintf(screen,"SNAP keyword %s %s \n",keywd,keyval);

      if (logfile) fprintf(logfile,"SNAP keyword %s %s \n",keywd,keyval);

    }



    if (strcmp(keywd,"rcutfac") == 0) {

      rcutfac = atof(keyval);

      rcutfacflag = 1;

    } else if (strcmp(keywd,"twojmax") == 0) {

      twojmax = atoi(keyval);

      twojmaxflag = 1;

    } else if (strcmp(keywd,"rfac0") == 0)

      rfac0 = atof(keyval);

    else if (strcmp(keywd,"rmin0") == 0)

      rmin0 = atof(keyval);

    else if (strcmp(keywd,"switchflag") == 0)

      switchflag = atoi(keyval);

    else if (strcmp(keywd,"bzeroflag") == 0)

      bzeroflag = atoi(keyval);

    else if (strcmp(keywd,"quadraticflag") == 0)

      quadraticflag = atoi(keyval);

    else if (strcmp(keywd,"alloyflag") == 0)

      alloyflag = atoi(keyval);

    else if (strcmp(keywd,"wselfallflag") == 0)

      wselfallflag = atoi(keyval);

    else if (strcmp(keywd,"chunksize") == 0)

      chunksize = atoi(keyval);

    else

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

  }



  bnormflag = alloyflag;



  if (rcutfacflag == 0 || twojmaxflag == 0)

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



}



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

   memory usage

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



double MLIAPDescriptorSNAP::memory_usage()

{

  double bytes = 0;



  int n = atom->ntypes+1;

  bytes += n*n*sizeof(int);      // setflag

  bytes += n*n*sizeof(double);   // cutsq

  bytes += n*sizeof(int);        // map

  bytes += beta_max*ncoeff*sizeof(double); // bispectrum

  bytes += beta_max*ncoeff*sizeof(double); // beta



  bytes += snaptr->memory_usage(); // SNA object



  return bytes;

}

src/MLIAP/mliap_descriptor_snap.h

0 → 100644
+67 −0
Original line number Diff line number Diff line 
/* -*- c++ -*- ----------------------------------------------------------

   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator

   http://lammps.sandia.gov, Sandia National Laboratories

   Steve Plimpton, sjplimp@sandia.gov



   Copyright (2003) Sandia Corporation.  Under the terms of Contract

   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains

   certain rights in this software.  This software is distributed under

   the GNU General Public License.



   See the README file in the top-level LAMMPS directory.

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



#ifndef MLIAP_DESCRIPTOR_SNAP_H

#define MLIAP_DESCRIPTOR_SNAP_H



#include "pair.h"



namespace LAMMPS_NS {



  //class MLIAPDescriptorSNAP : protected Pointers  {

class MLIAPDescriptorSNAP : public Pair  {

public:

  MLIAPDescriptorSNAP(class LAMMPS *);

  ~MLIAPDescriptorSNAP();

  virtual void compute(int, int);

  void settings(int, char **);

  virtual void coeff(int, char **);

  virtual void init_style();

  virtual double init_one(int, int);

  virtual double memory_usage();



  double rcutfac, quadraticflag; // declared public to workaround gcc 4.9

  int ncoeff;                    //  compiler bug, manifest in KOKKOS package



protected:

  int ncoeffq, ncoeffall;

  class SNA* snaptr;

  virtual void allocate();

  void read_files(char *, char *);

  inline int equal(double* x,double* y);

  inline double dist2(double* x,double* y);



  void compute_beta();

  void compute_bispectrum();



  double rcutmax;               // max cutoff for all elements

  int nelements;                // # of unique elements

  char **elements;              // names of unique elements

  double *radelem;              // element radii

  double *wjelem;               // elements weights

  double **coeffelem;           // element bispectrum coefficients

  double** beta;                // betas for all atoms in list

  double** bispectrum;          // bispectrum components for all atoms in list

  int *map;                     // mapping from atom types to elements

  int twojmax, switchflag, bzeroflag, bnormflag;

  int alloyflag, wselfallflag;

  int chunksize;

  double rfac0, rmin0, wj1, wj2;

  int rcutfacflag, twojmaxflag; // flags for required parameters

  int beta_max;                 // length of beta

};



}



#endif

src/MLIAP/mliap_model_snap.cpp

0 → 100644
+426 −0

File added.

Preview size limit exceeded, changes collapsed.

src/MLIAP/mliap_model_snap.h

0 → 100644
+59 −0

File added.

Preview size limit exceeded, changes collapsed.

src/MLIAP/pair_mliap_snap.cpp

0 → 100644
+702 −0

File added.

Preview size limit exceeded, changes collapsed.

CRA Git | Maintained and supported by SUSTech CRA and CCSE