diff --git a/src/USER-MISC/README b/src/USER-MISC/README index 79b6e4373..ba19f77ac 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -1,53 +1,52 @@ The files in this package are a potpourri of (mostly) unrelated features contributed to LAMMPS by users. Each feature is a single pair of files (*.cpp and *.h). More information about each feature can be found by reading its doc page in the LAMMPS doc directory. The doc page which lists all LAMMPS input script commands is as follows: doc/Section_commands.html, subsection 3.5 User-contributed features are listed at the bottom of the fix, compute, pair, etc sections. The list of features and author of each is given below. You should contact the author directly if you have specific questions about the feature or its coding. ------------------------------------------------------------ angle_style cosine/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 angle_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 angle_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 angle_style fourier/simple, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 angle_style dipole, Mario Orsi, orsimario at gmail.com, 10 Jan 12 angle_style quartic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007 compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013 compute temp/rotate, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 dihedral_style cosine/shift/exp, Carsten Svaneborg, science at zqex.dk, 8 Aug 11 dihedral_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 dihedral_style nharmonic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 dihedral_style quadratic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 fix imd, Axel Kohlmeyer, akohlmey at gmail.com, 9 Nov 2009 fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008 improper_style cossq, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12 improper_style fourier, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12 improper_style ring, Georgios Vogiatzis, gvog at chemeng.ntua.gr, 25 May 12 pair_style coul/diel, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11 pair_style dipole/sf, Mario Orsi, orsimario at gmail.com, 8 Aug 11 pair_style edip, Luca Ferraro, luca.ferraro at caspur.it, 15 Sep 11 pair_style eam/cd, Alexander Stukowski, stukowski at mm.tu-darmstadt.de, 7 Nov 09 pair_style gauss/cut, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11 pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13 pair_style lj/sf, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12 pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12 -pair_style nb3b/harmonic, Todd R. Zeitler, todd.zeitler at gmail.com, 6 Apr 13 pair_style tersoff/table, Luca Ferraro, luca.ferraro@caspur.it, 1 Dec 11 diff --git a/src/USER-MISC/pair_nb3b_harmonic.cpp b/src/USER-MISC/pair_nb3b_harmonic.cpp deleted file mode 100644 index 1130479d3..000000000 --- a/src/USER-MISC/pair_nb3b_harmonic.cpp +++ /dev/null @@ -1,532 +0,0 @@ -/* ---------------------------------------------------------------------- - 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: Todd R. Zeitler (SNL) - (based on Stillinger-Weber pair style) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "pair_nb3b_harmonic.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 -#define SMALL 0.001 -#define PI 3.141592653589793238462643383279 - -/* ---------------------------------------------------------------------- */ - -PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 0; - one_coeff = 1; - manybody_flag = 1; - - nelements = 0; - elements = NULL; - nparams = maxparam = 0; - params = NULL; - elem2param = NULL; -} - -/* ---------------------------------------------------------------------- - check if allocated, since class can be destructed when incomplete -------------------------------------------------------------------------- */ - -PairNb3bHarmonic::~PairNb3bHarmonic() -{ - if (elements) - for (int i = 0; i < nelements; i++) delete [] elements[i]; - delete [] elements; - memory->destroy(params); - memory->destroy(elem2param); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - delete [] map; - } -} - -/* ---------------------------------------------------------------------- */ - -void PairNb3bHarmonic::compute(int eflag, int vflag) -{ - int i,j,k,ii,jj,kk,inum,jnum,jnumm1,itag,jtag; - int itype,jtype,ktype,ijparam,ikparam,ijkparam; - 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; - int *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]; - - // 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]; - - 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; - } - - 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].cutsq) continue; - - } - - 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].cutsq) 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].cutsq) 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 PairNb3bHarmonic::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 PairNb3bHarmonic::settings(int narg, char **arg) -{ - if (narg != 0) error->all(FLERR,"Illegal pair_style command"); -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairNb3bHarmonic::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(); - - // 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 PairNb3bHarmonic::init_style() -{ - if (atom->tag_enable == 0) - error->all(FLERR,"Pair style nb3b/harmonic requires atom IDs"); - if (force->newton_pair == 0) - error->all(FLERR,"Pair style nb3b/harmonic 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 PairNb3bHarmonic::init_one(int i, int j) -{ - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); - - return cutmax; -} - -/* ---------------------------------------------------------------------- */ - -void PairNb3bHarmonic::read_file(char *file) -{ - int params_per_line = 6; - 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 = open_potential(file); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open nb3b/harmonic 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 nb3b/harmonic 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].k_theta = atof(words[3]); - params[nparams].theta0 = atof(words[4]); - params[nparams].cutoff = atof(words[5]); - - if (params[nparams].k_theta < 0.0 || params[nparams].theta0 < 0.0 || - params[nparams].cutoff < 0.0) - error->all(FLERR,"Illegal nb3b/harmonic parameter"); - - nparams++; - } - - delete [] words; -} - -/* ---------------------------------------------------------------------- */ - -void PairNb3bHarmonic::setup() -{ - int i,j,k,m,n; - double rtmp; - - // 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 - // (cut = a*sigma) - - for (m = 0; m < nparams; m++) { - - params[m].cut = params[m].cutoff; - params[m].cutsq = params[m].cut * params[m].cut; - - params[m].theta0 = params[m].theta0 / 180.0 * PI; - - } - - // set cutmax to max of all params - - cutmax = 0.0; - for (m = 0; m < nparams; m++) { - rtmp = sqrt(params[m].cutsq); - if (rtmp > cutmax) cutmax = rtmp; - } -} - -/* ---------------------------------------------------------------------- */ - - -void PairNb3bHarmonic::threebody(Param *paramij, Param *paramik, - Param *paramijk, - double rsq1, double rsq2, - double *delr1, double *delr2, - double *fj, double *fk, int eflag, double &eng) -{ - double dtheta,tk; - double r1,r2,c,s,a,a11,a12,a22; - - // angle (cos and sin) - - r1 = sqrt(rsq1); - r2 = sqrt(rsq2); - - c = delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]; - c /= r1*r2; - - if (c > 1.0) c = 1.0; - if (c < -1.0) c = -1.0; - - s = sqrt(1.0 - c*c); - if (s < SMALL) s = SMALL; - s = 1.0/s; - - // force & energy - - dtheta = acos(c) - paramijk->theta0; - tk = paramijk->k_theta * dtheta; - - if (eflag) eng = tk*dtheta; - - a = -2.0 * tk * s; - a11 = a*c / rsq1; - a12 = -a / (r1*r2); - a22 = a*c / rsq2; - - fj[0] = a11*delr1[0] + a12*delr2[0]; - fj[1] = a11*delr1[1] + a12*delr2[1]; - fj[2] = a11*delr1[2] + a12*delr2[2]; - fk[0] = a22*delr2[0] + a12*delr1[0]; - fk[1] = a22*delr2[1] + a12*delr1[1]; - fk[2] = a22*delr2[2] + a12*delr1[2]; -} - diff --git a/src/USER-MISC/pair_nb3b_harmonic.h b/src/USER-MISC/pair_nb3b_harmonic.h deleted file mode 100644 index 5d126ae09..000000000 --- a/src/USER-MISC/pair_nb3b_harmonic.h +++ /dev/null @@ -1,64 +0,0 @@ -/* --*- 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. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(nb3b/harmonic,PairNb3bHarmonic) - -#else - -#ifndef LMP_PAIR_NB3B_HARMONIC_H -#define LMP_PAIR_NB3B_HARMONIC_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairNb3bHarmonic : public Pair { - public: - PairNb3bHarmonic(class LAMMPS *); - virtual ~PairNb3bHarmonic(); - virtual void compute(int, int); - void settings(int, char **); - void coeff(int, char **); - double init_one(int, int); - void init_style(); - - protected: - struct Param { - double k_theta, theta0, cutoff; - double cut,cutsq; - int ielement,jelement,kelement; - }; - - 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 - int nparams; // # of stored parameter sets - int maxparam; // max # of parameter sets - Param *params; // parameter set for an I-J-K interaction - - void allocate(); - void read_file(char *); - void setup(); - void twobody(Param *, double, double &, int, double &); - void threebody(Param *, Param *, Param *, double, double, double *, double *, - double *, double *, int, double &); -}; - -} - -#endif -#endif