From a8ef38ce7d14a237aeb651325f75f4b9f0422a38 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 29 Oct 2013 15:16:11 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10881 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/USER-MISC/README | 1 - src/USER-MISC/pair_nb3b_harmonic.cpp | 532 --------------------------- src/USER-MISC/pair_nb3b_harmonic.h | 64 ---- 3 files changed, 597 deletions(-) delete mode 100644 src/USER-MISC/pair_nb3b_harmonic.cpp delete mode 100644 src/USER-MISC/pair_nb3b_harmonic.h diff --git a/src/USER-MISC/README b/src/USER-MISC/README index 79b6e4373d..ba19f77ac0 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -49,5 +49,4 @@ 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 1130479d3f..0000000000 --- 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 5d126ae09a..0000000000 --- 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