diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index fa2be92b7..5c5e4804f 100755 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -1,729 +1,741 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Contributing author: Yongnan Xiong (HNU), xyn@hnu.edu.cn Aidan Thompson (SNL) Anders Hafreager (UiO), andershaf@gmail.com ------------------------------------------------------------------------- */ #include #include #include #include #include "pair_vashishta.h" #include "atom.h" #include "neighbor.h" #include "neigh_request.h" #include "force.h" #include "comm.h" #include "memory.h" #include "neighbor.h" #include "neigh_list.h" #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; #define MAXLINE 1024 #define DELTA 4 /* ---------------------------------------------------------------------- */ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; restartinfo = 0; one_coeff = 1; manybody_flag = 1; nelements = 0; elements = NULL; nparams = maxparam = 0; params = NULL; elem2param = NULL; map = NULL; neigh3BodyMax = 0; neigh3BodyCount = NULL; neigh3Body = NULL; - useTable = true; - tableSize = 65536; + useTable = false; + nTablebits = 12; forceTable = NULL; potentialTable = NULL; } /* ---------------------------------------------------------------------- check if allocated, since class can be destructed when incomplete ------------------------------------------------------------------------- */ PairVashishta::~PairVashishta() { if (elements) for (int i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; memory->destroy(params); memory->destroy(elem2param); memory->destroy(forceTable); memory->destroy(potentialTable); memory->destroy(neigh3BodyCount); memory->destroy(neigh3Body); if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); delete [] map; } } void PairVashishta::modify_params(int narg, char **arg) { - if(narg < 2 || narg > 3) // We accept table yes/no [tableSize] - error->all(FLERR,"Illegal pair_modify command"); - if(strcmp(arg[0], "table") != 0) - error->all(FLERR,"Illegal pair_modify command. Expected 'table yes/no [tableSize]'"); - if(strcmp(arg[1], "yes") == 0) - useTable = true; - else if(strcmp(arg[1], "no") == 0) - useTable = false; - else - error->all(FLERR,"Illegal pair_modify command. Expected 'table yes/no [tableSize]'"); - - if(narg == 3) { - tableSize = atoi(arg[2]); - if(tableSize <= 0) - error->all(FLERR,"Illegal pair_modify command. Expected 'table yes/no [tableSize]', but tableSize has value 0."); - - createTable(); + if (narg == 0) error->all(FLERR,"Illegal pair_modify command"); + + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg],"table") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); + nTablebits = force->inumeric(FLERR,arg[iarg+1]); + if (nTablebits > sizeof(float)*CHAR_BIT) + error->all(FLERR,"Too many total bits for bitmapped lookup table"); + if(nTablebits != 0) + useTable = true; + iarg += 2; + } else if (strcmp(arg[iarg],"tabinner") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); + tabinner = force->numeric(FLERR,arg[iarg+1]); + iarg += 2; + } } + + createTable(); } void PairVashishta::validateNeigh3Body() { if (atom->nlocal > neigh3BodyMax) { neigh3BodyMax = atom->nmax; memory->destroy(neigh3BodyCount); memory->destroy(neigh3Body); memory->create(neigh3BodyCount,neigh3BodyMax, "pair:vashishta:neigh3BodyCount"); memory->create(neigh3Body,neigh3BodyMax, 1000, "pair:vashishta:neigh3Body"); // TODO: pick this number more wisely? } } void PairVashishta::createTable() { int ntypes = atom->ntypes+1; - deltaR2 = cutmax*cutmax / (tableSize-1); + tabinnersq = tabinner*tabinner; + + int ntable = 1; + for (int i = 0; i < nTablebits; i++) ntable *= 2; + + deltaR2 = (cutmax*cutmax - tabinnersq) / (ntable-1); oneOverDeltaR2 = 1.0/deltaR2; memory->destroy(forceTable); memory->destroy(potentialTable); - memory->create(forceTable,ntypes,ntypes,tableSize+1,"pair:vashishta:forceTable"); - memory->create(potentialTable,ntypes,ntypes,tableSize+1,"pair:vashishta:potentialTable"); + memory->create(forceTable,ntypes,ntypes,ntable+1,"pair:vashishta:forceTable"); + memory->create(potentialTable,ntypes,ntypes,ntable+1,"pair:vashishta:potentialTable"); for (int ii = 0; ii < ntypes; ii++) { int i = map[ii]; for (int jj = 0; jj < ntypes; jj++) { - int j = map[jj]; - if (i < 0 || j < 0 || ii == 0 || jj == 0) { - for(int tableIndex=0; tableIndex<=tableSize; tableIndex++) { - forceTable[ii][jj][tableIndex] = 0; - potentialTable[ii][jj][tableIndex] = 0; - } - } else { - int ijparam = elem2param[i][j][j]; - for(int tableIndex=0; tableIndex<=tableSize; tableIndex++) { - double rsq = tableIndex*deltaR2; - double fpair, eng; - twobody(¶ms[ijparam], rsq, fpair, 1, eng, false /* Don't ask for tabulated since we are generating it now*/ ); - forceTable[ii][jj][tableIndex] = fpair; - potentialTable[ii][jj][tableIndex] = eng; - } + int j = map[jj]; + if (i < 0 || j < 0 || ii == 0 || jj == 0) { + for(int tableIndex=0; tableIndex<=ntable; tableIndex++) { + forceTable[ii][jj][tableIndex] = 0; + potentialTable[ii][jj][tableIndex] = 0; } + } else { + int ijparam = elem2param[i][j][j]; + for(int tableIndex=0; tableIndex<=ntable; tableIndex++) { + double rsq = tabinnersq + tableIndex*deltaR2; + double fpair, eng; + twobody(¶ms[ijparam], rsq, fpair, 1, eng, false /* Don't ask for tabulated since we are generating it now*/ ); + forceTable[ii][jj][tableIndex] = fpair; + potentialTable[ii][jj][tableIndex] = eng; + } + } } } } /* ---------------------------------------------------------------------- */ void PairVashishta::compute(int eflag, int vflag) { validateNeigh3Body(); int i,j,k,ii,jj,kk,inum,jnum,jnumm1; int itype,jtype,ktype,ijparam,ikparam,ijkparam; tagint itag,jtag; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,rsq1,rsq2; double delr1[3],delr2[3],fj[3],fk[3]; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; tagint *tag = atom->tag; int *type = atom->type; int nlocal = atom->nlocal; int newton_pair = force->newton_pair; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; // loop over full neighbor list of my atoms for (ii = 0; ii < inum; ii++) { i = ilist[ii]; itag = tag[i]; itype = map[type[i]]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; neigh3BodyCount[i] = 0; // Reset the 3-body neighbor list // two-body interactions, skip half of them jlist = firstneigh[i]; jnum = numneigh[i]; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; jtag = tag[j]; jtype = map[type[j]]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; ijparam = elem2param[itype][jtype][jtype]; if (rsq <= params[ijparam].cutsq2) { neigh3Body[i][neigh3BodyCount[i]] = j; neigh3BodyCount[i]++; } if (rsq > params[ijparam].cutsq) continue; if (itag > jtag) { if ((itag+jtag) % 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) % 2 == 1) continue; } else { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp && x[j][1] < ytmp) continue; if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue; } twobody(¶ms[ijparam],rsq,fpair,eflag,evdwl, useTable); f[i][0] += delx*fpair; f[i][1] += dely*fpair; f[i][2] += delz*fpair; f[j][0] -= delx*fpair; f[j][1] -= dely*fpair; f[j][2] -= delz*fpair; if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,0.0,fpair,delx,dely,delz); } jlist = neigh3Body[i]; jnum = neigh3BodyCount[i]; jnumm1 = jnum - 1; for (jj = 0; jj < jnumm1; jj++) { j = jlist[jj]; j &= NEIGHMASK; jtype = map[type[j]]; ijparam = elem2param[itype][jtype][jtype]; delr1[0] = x[j][0] - xtmp; delr1[1] = x[j][1] - ytmp; delr1[2] = x[j][2] - ztmp; rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2]; if (rsq1 >= params[ijparam].cutsq2) continue; for (kk = jj+1; kk < jnum; kk++) { k = jlist[kk]; k &= NEIGHMASK; ktype = map[type[k]]; ikparam = elem2param[itype][ktype][ktype]; ijkparam = elem2param[itype][jtype][ktype]; delr2[0] = x[k][0] - xtmp; delr2[1] = x[k][1] - ytmp; delr2[2] = x[k][2] - ztmp; rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2]; if (rsq2 >= params[ikparam].cutsq2) continue; threebody(¶ms[ijparam],¶ms[ikparam],¶ms[ijkparam], rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl); f[i][0] -= fj[0] + fk[0]; f[i][1] -= fj[1] + fk[1]; f[i][2] -= fj[2] + fk[2]; f[j][0] += fj[0]; f[j][1] += fj[1]; f[j][2] += fj[2]; f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2]; if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2); } } } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- */ void PairVashishta::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); map = new int[n+1]; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairVashishta::settings(int narg, char **arg) { if (narg != 0) error->all(FLERR,"Illegal pair_style command"); } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairVashishta::coeff(int narg, char **arg) { int i,j,n; if (!allocated) allocate(); if (narg != 3 + atom->ntypes) error->all(FLERR,"Incorrect args for pair coefficients"); // insure I,J args are * * if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->all(FLERR,"Incorrect args for pair coefficients"); // read args that map atom types to elements in potential file // map[i] = which element the Ith atom type is, -1 if NULL // nelements = # of unique elements // elements = list of element names if (elements) { for (i = 0; i < nelements; i++) delete [] elements[i]; delete [] elements; } elements = new char*[atom->ntypes]; for (i = 0; i < atom->ntypes; i++) elements[i] = NULL; nelements = 0; for (i = 3; i < narg; i++) { if (strcmp(arg[i],"NULL") == 0) { map[i-2] = -1; continue; } for (j = 0; j < nelements; j++) if (strcmp(arg[i],elements[j]) == 0) break; map[i-2] = j; if (j == nelements) { n = strlen(arg[i]) + 1; elements[j] = new char[n]; strcpy(elements[j],arg[i]); nelements++; } } // read potential file and initialize potential parameters read_file(arg[2]); setup_params(); // clear setflag since coeff() called once with I,J = * * 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"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairVashishta::init_style() { if (atom->tag_enable == 0) error->all(FLERR,"Pair style Vashishta requires atom IDs"); if (force->newton_pair == 0) error->all(FLERR,"Pair style Vashishta requires newton pair on"); // need a full neighbor list int irequest = neighbor->request(this); neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->full = 1; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairVashishta::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); return cutmax; } /* ---------------------------------------------------------------------- */ void PairVashishta::read_file(char *file) { int params_per_line = 17; char **words = new char*[params_per_line+1]; memory->sfree(params); params = NULL; nparams = maxparam = 0; // open file on proc 0 FILE *fp; if (comm->me == 0) { fp = force->open_potential(file); if (fp == NULL) { char str[128]; sprintf(str,"Cannot open Vashishta potential file %s",file); error->one(FLERR,str); } } // read each set of params from potential file // one set of params can span multiple lines // store params if all 3 element tags are in element list int n,nwords,ielement,jelement,kelement; char line[MAXLINE],*ptr; int eof = 0; while (1) { if (comm->me == 0) { ptr = fgets(line,MAXLINE,fp); if (ptr == NULL) { eof = 1; fclose(fp); } 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; // concatenate additional lines until have params_per_line words while (nwords < params_per_line) { n = strlen(line); if (comm->me == 0) { ptr = fgets(&line[n],MAXLINE-n,fp); if (ptr == NULL) { eof = 1; fclose(fp); } 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); if ((ptr = strchr(line,'#'))) *ptr = '\0'; nwords = atom->count_words(line); } if (nwords != params_per_line) error->all(FLERR,"Incorrect format in Vashishta potential file"); // words = ptrs to all words in line nwords = 0; words[nwords++] = strtok(line," \t\n\r\f"); while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue; // ielement,jelement,kelement = 1st args // if all 3 args are in element list, then parse this line // else skip to next entry in file for (ielement = 0; ielement < nelements; ielement++) if (strcmp(words[0],elements[ielement]) == 0) break; if (ielement == nelements) continue; for (jelement = 0; jelement < nelements; jelement++) if (strcmp(words[1],elements[jelement]) == 0) break; if (jelement == nelements) continue; for (kelement = 0; kelement < nelements; kelement++) if (strcmp(words[2],elements[kelement]) == 0) break; if (kelement == nelements) continue; // load up parameter settings and error check their values if (nparams == maxparam) { maxparam += DELTA; params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params"); } params[nparams].ielement = ielement; params[nparams].jelement = jelement; params[nparams].kelement = kelement; params[nparams].bigh = atof(words[3]); params[nparams].eta = atof(words[4]); params[nparams].zi = atof(words[5]); params[nparams].zj = atof(words[6]); params[nparams].lambda1 = atof(words[7]); params[nparams].bigd = atof(words[8]); params[nparams].lambda4 = atof(words[9]); params[nparams].bigw = atof(words[10]); params[nparams].cut = atof(words[11]); params[nparams].bigb = atof(words[12]); params[nparams].gamma = atof(words[13]); params[nparams].r0 = atof(words[14]); params[nparams].bigc = atof(words[15]); params[nparams].costheta = atof(words[16]); if (params[nparams].bigb < 0.0 || params[nparams].gamma < 0.0 || params[nparams].r0 < 0.0 || params[nparams].bigc < 0.0 || params[nparams].bigh < 0.0 || params[nparams].eta < 0.0 || params[nparams].lambda1 < 0.0 || params[nparams].bigd < 0.0 || params[nparams].lambda4 < 0.0 || params[nparams].bigw < 0.0 || params[nparams].cut < 0.0) error->all(FLERR,"Illegal Vashishta parameter"); nparams++; } delete [] words; } /* ---------------------------------------------------------------------- */ void PairVashishta::setup_params() { int i,j,k,m,n; // set elem2param for all triplet combinations // must be a single exact match to lines read from file // do not allow for ACB in place of ABC memory->destroy(elem2param); memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param"); for (i = 0; i < nelements; i++) for (j = 0; j < nelements; j++) for (k = 0; k < nelements; k++) { n = -1; for (m = 0; m < nparams; m++) { if (i == params[m].ielement && j == params[m].jelement && k == params[m].kelement) { if (n >= 0) error->all(FLERR,"Potential file has duplicate entry"); n = m; } } if (n < 0) error->all(FLERR,"Potential file is missing an entry"); elem2param[i][j][k] = n; } // compute parameter values derived from inputs // set cutsq using shortcut to reduce neighbor list for accelerated // calculations. cut must remain unchanged as it is a potential parameter double tmp_par; for (m = 0; m < nparams; m++) { params[m].cutsq = params[m].cut * params[m].cut; params[m].cutsq2 = params[m].r0 * params[m].r0; tmp_par = params[m].lambda1; params[m].lam1inv = (tmp_par == 0.0) ? 0.0 : 1.0/tmp_par; tmp_par = params[m].lambda4; params[m].lam4inv = (tmp_par == 0.0) ? 0.0 : 1.0/tmp_par; params[m].zizj = params[m].zi*params[m].zj * force->qqr2e; // note that bigd does not have 1/2 factor params[m].mbigd = params[m].bigd; params[m].heta = params[m].bigh*params[m].eta; params[m].big2b = 2.0*params[m].bigb; params[m].big6w = 6.0*params[m].bigw; tmp_par = params[m].cut; params[m].rcinv = (tmp_par == 0.0) ? 0.0 : 1.0/tmp_par; params[m].rc2inv = params[m].rcinv*params[m].rcinv; params[m].rc4inv = params[m].rc2inv*params[m].rc2inv; params[m].rc6inv = params[m].rc2inv*params[m].rc4inv; params[m].rceta = pow(params[m].rcinv,params[m].eta); params[m].lam1rc = params[m].cut*params[m].lam1inv; params[m].lam4rc = params[m].cut*params[m].lam4inv; params[m].vrcc2 = params[m].zizj*params[m].rcinv * exp(-params[m].lam1rc); params[m].vrcc3 = params[m].mbigd*params[m].rc4inv * exp(-params[m].lam4rc); params[m].vrc = params[m].bigh*params[m].rceta + params[m].vrcc2 - params[m].vrcc3 - params[m].bigw*params[m].rc6inv; params[m].dvrc = params[m].vrcc3 * (4.0*params[m].rcinv+params[m].lam4inv) + params[m].big6w * params[m].rc6inv * params[m].rcinv - params[m].heta * params[m].rceta*params[m].rcinv - params[m].vrcc2 * (params[m].rcinv+params[m].lam1inv); params[m].c0 = params[m].cut*params[m].dvrc - params[m].vrc; } // set cutmax to max of all params cutmax = 0.0; for (m = 0; m < nparams; m++) { if (params[m].cut > cutmax) cutmax = params[m].cut; if (params[m].r0 > cutmax) cutmax = params[m].r0; } createTable(); } /* ---------------------------------------------------------------------- */ void PairVashishta::twobody(Param *param, double rsq, double &fforce, int eflag, double &eng, bool tabulated) { if(tabulated) { - int potentialTableIndex = rsq*oneOverDeltaR2; - double fractionalRest = rsq*oneOverDeltaR2 - potentialTableIndex; // double - int will only keep the 0.xxxx part + if (rsq < tabinnersq) { + sprintf(estr,"Pair distance < table inner cutoff: " + "ijtype %d %d dist %g",param->ielement+1,param->jelement+1,sqrt(rsq)); + error->one(FLERR,estr); + } + + int tableIndex = (rsq - tabinnersq)*oneOverDeltaR2; + double fraction = (rsq - tabinnersq)*oneOverDeltaR2 - tableIndex; // double - int will only keep the 0.xxxx part - double force0 = forceTable[param->ielement+1][param->jelement+1][potentialTableIndex]; - double force1 = forceTable[param->ielement+1][param->jelement+1][potentialTableIndex+1]; - fforce = (1.0 - fractionalRest)*force0 + fractionalRest*force1; // force is linearly interpolated between the two values + double force0 = forceTable[param->ielement+1][param->jelement+1][tableIndex]; + double force1 = forceTable[param->ielement+1][param->jelement+1][tableIndex+1]; + fforce = (1.0 - fraction)*force0 + fraction*force1; // force is linearly interpolated between the two values if(evflag) { - double energy0 = potentialTable[param->ielement+1][param->jelement+1][potentialTableIndex]; - double energy1 = potentialTable[param->ielement+1][param->jelement+1][potentialTableIndex+1]; - eng = (1.0 - fractionalRest)*energy0 + fractionalRest*energy1; + double energy0 = potentialTable[param->ielement+1][param->jelement+1][tableIndex]; + double energy1 = potentialTable[param->ielement+1][param->jelement+1][tableIndex+1]; + eng = (1.0 - fraction)*energy0 + fraction*energy1; } } else { double r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r,vc2,vc3; r = sqrt(rsq); rinvsq = 1.0/rsq; r4inv = rinvsq*rinvsq; r6inv = rinvsq*r4inv; reta = pow(r,-param->eta); lam1r = r*param->lam1inv; lam4r = r*param->lam4inv; vc2 = param->zizj * exp(-lam1r)/r; vc3 = param->mbigd * r4inv*exp(-lam4r); fforce = (param->dvrc*r - (4.0*vc3 + lam4r*vc3+param->big6w*r6inv - param->heta*reta - vc2 - lam1r*vc2) ) * rinvsq; if (eflag) eng = param->bigh*reta + vc2 - vc3 - param->bigw*r6inv - r*param->dvrc + param->c0; } } /* ---------------------------------------------------------------------- */ void PairVashishta::threebody(Param *paramij, Param *paramik, Param *paramijk, double rsq1, double rsq2, double *delr1, double *delr2, double *fj, double *fk, int eflag, double &eng) { double r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1; double r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2; double rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2,pcsinv,pcsinvsq,pcs; double facang,facang12,csfacang,csfac1,csfac2; r1 = sqrt(rsq1); rinvsq1 = 1.0/rsq1; rainv1 = 1.0/(r1 - paramij->r0); gsrainv1 = paramij->gamma * rainv1; gsrainvsq1 = gsrainv1*rainv1/r1; expgsrainv1 = exp(gsrainv1); r2 = sqrt(rsq2); rinvsq2 = 1.0/rsq2; rainv2 = 1.0/(r2 - paramik->r0); gsrainv2 = paramik->gamma * rainv2; gsrainvsq2 = gsrainv2*rainv2/r2; expgsrainv2 = exp(gsrainv2); rinv12 = 1.0/(r1*r2); cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12; delcs = cs - paramijk->costheta; delcssq = delcs*delcs; pcsinv = paramijk->bigc*delcssq + 1.0; pcsinvsq = pcsinv*pcsinv; pcs = delcssq/pcsinv; facexp = expgsrainv1*expgsrainv2; facrad = paramijk->bigb * facexp * pcs; frad1 = facrad*gsrainvsq1; frad2 = facrad*gsrainvsq2; facang = paramijk->big2b * facexp * delcs/pcsinvsq; facang12 = rinv12*facang; csfacang = cs*facang; csfac1 = rinvsq1*csfacang; fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12; fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12; fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12; csfac2 = rinvsq2*csfacang; fk[0] = delr2[0]*(frad2+csfac2)-delr1[0]*facang12; fk[1] = delr2[1]*(frad2+csfac2)-delr1[1]*facang12; fk[2] = delr2[2]*(frad2+csfac2)-delr1[2]*facang12; if (eflag) eng = facrad; } diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h index 47a36e39a..47c4ab9b7 100755 --- a/src/MANYBODY/pair_vashishta.h +++ b/src/MANYBODY/pair_vashishta.h @@ -1,132 +1,134 @@ /* ---------------------------------------------------------------------- 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. ------------------------------------------------------------------------- */ #ifdef PAIR_CLASS PairStyle(vashishta,PairVashishta) #else #ifndef LMP_PAIR_Vashishta_H #define LMP_PAIR_Vashishta_H #include "pair.h" namespace LAMMPS_NS { class PairVashishta : public Pair { public: PairVashishta(class LAMMPS *); virtual ~PairVashishta(); virtual void modify_params(int, char **); virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); virtual double init_one(int, int); virtual void init_style(); protected: struct Param { double bigb,gamma,r0,bigc,costheta; double bigh,eta,zi,zj; double lambda1,bigd,mbigd,lambda4,bigw,cut; double lam1inv,lam4inv,zizj,heta,big2b,big6w; double rcinv,rc2inv,rc4inv,rc6inv,rceta; double cutsq2,cutsq; double lam1rc,lam4rc,vrcc2,vrcc3,vrc,dvrc,c0; int ielement,jelement,kelement; }; + bool useTable; - int tableSize; + int nTablebits; double deltaR2; double oneOverDeltaR2; double ***forceTable; // Table containing the forces double ***potentialTable; // Table containing the potential energies int neigh3BodyMax; // Max size of short neighborlist int *neigh3BodyCount; // Number of neighbors in short range 3 particle forces neighbor list int **neigh3Body; // Neighborlist for short range 3 particle forces double cutmax; // max cutoff for all elements int nelements; // # of unique elements char **elements; // names of unique elements int ***elem2param; // mapping from element triplets to parameters int *map; // mapping from atom types to elements + char estr[128]; // Char array to store error message with sprintf (i.e. twobody table) int nparams; // # of stored parameter sets int maxparam; // max # of parameter sets Param *params; // parameter set for an I-J-K interaction virtual void allocate(); void read_file(char *); void setup_params(); void createTable(); void validateNeigh3Body(); - void twobody(Param *, double, double &, int, double &, bool); + void twobody(Param *, double, double &, int, double &, bool tabulated = false); void threebody(Param *, Param *, Param *, double, double, double *, double *, double *, double *, int, double &); }; } #endif #endif /* ERROR/WARNING messages: E: Illegal ... command Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. E: Incorrect args for pair coefficients Self-explanatory. Check the input script or data file. E: Pair style Vashishta requires atom IDs This is a requirement to use the Vashishta potential. E: Pair style Vashishta requires newton pair on See the newton command. This is a restriction to use the Vashishta potential. E: All pair coeffs are not set All pair coefficients must be set in the data file or by the pair_coeff command before running a simulation. E: Cannot open Vashishta potential file %s The specified Vashishta potential file cannot be opened. Check that the path and name are correct. E: Incorrect format in Vashishta potential file Incorrect number of words per line in the potential file. E: Illegal Vashishta parameter One or more of the coefficients defined in the potential file is invalid. E: Potential file has duplicate entry The potential file has more than one entry for the same element. E: Potential file is missing an entry The potential file does not have a needed entry. */