diff --git a/src/MANYBODY/pair_polymorphic.cpp b/src/MANYBODY/pair_polymorphic.cpp index 0d3587204..db99a86b3 100755 --- a/src/MANYBODY/pair_polymorphic.cpp +++ b/src/MANYBODY/pair_polymorphic.cpp @@ -1,796 +1,978 @@ -/* ---------------------------------------------------------------------- - 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 authors: Xiaowang Zhou, Reese Jones (SNL) - Based on pair_tersoff by Aidan Thompson (SNL) -------------------------------------------------------------------------- */ - -#include -#include -#include -#include -#include "pair_polymorphic.h" -#include "atom.h" -#include "neighbor.h" -#include "neigh_list.h" -#include "neigh_request.h" -#include "force.h" -#include "comm.h" -#include "memory.h" -#include "error.h" - -#include "math_const.h" - -using namespace LAMMPS_NS; -using namespace MathConst; - -#define MAXLINE 1024 -#define DELTA 4 - -/* ====================================================================== */ - -PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp) -{ - single_enable = 0; - one_coeff = 1; - - nelements = 0; - elements = NULL; - pairParameters = NULL; - tripletParameters = NULL; - elem2param = NULL; - elem3param = NULL; - type_map = NULL; -} - -/* ---------------------------------------------------------------------- - check if allocated, since class can be destructed when incomplete -------------------------------------------------------------------------- */ - -PairPolymorphic::~PairPolymorphic() -{ - if (elements) - for (int i = 0; i < nelements; i++) delete [] elements[i]; - delete [] elements; - memory->destroy(pairParameters); - memory->destroy(tripletParameters); - memory->destroy(elem2param); - memory->destroy(elem3param); - - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - delete [] type_map; - } -} - -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::compute(int eflag, int vflag) -{ - tagint itag,jtag; - int i,j,k,ii,jj,kk,inum,jnum; - int iel,jel,kel,iparam_ij,iparam_ik,iparam_ijk; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,rsq1,rsq2,r0,r1,r2; - double delr1[3],delr2[3],fi[3],fj[3],fk[3]; - double zeta_ij,prefactor,wfac,pfac,gfac,fa,fa_d,bij,bij_d; - double costheta; - int *ilist,*jlist,*numneigh,**firstneigh; - - evdwl = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = vflag_atom = 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]; - iel = type_map[type[i]]; - if (iel < 0) continue; - 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] < x[i][2]) 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; - } - - - jel = type_map[type[j]]; - if (jel < 0) continue; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - - iparam_ij = elem2param[iel][jel]; - PairParameters & p = pairParameters[iparam_ij]; - if (rsq > p.cutsq) continue; - r0 = sqrt(rsq); - if (eflag) evdwl = (p.U)->value(r0); - fpair = (p.U)->derivative(r0); - fpair = -fpair/r0; - - 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); - } - - if (eta) { - - iparam_ij = elem2param[iel][iel]; - PairParameters & p = pairParameters[iparam_ij]; - - // accumulate bondorder zeta for each i-j interaction via loop over k - - zeta_ij = 0.0; - - for (kk = 0; kk < jnum; kk++) { - k = jlist[kk]; - k &= NEIGHMASK; - kel = type_map[type[k]]; - if (kel < 0) continue; - - 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]; - - iparam_ik = elem2param[kel][kel]; - PairParameters & q = pairParameters[iparam_ik]; - if (rsq2 > q.cutsq) continue; - r2 = sqrt(rsq2); - - wfac = (q.W)->value(r2); - - zeta_ij += wfac; - } - - // pairwise force due to zeta - - bij = (p.F)->value(zeta_ij); - bij_d = (p.F)->derivative(zeta_ij); - - prefactor = 0.5* bij_d; - if (eflag) evdwl = -0.5*bij; - - if (evflag) ev_tally(i,i,nlocal,newton_pair, - evdwl,0.0,0.0,-delr1[0],-delr1[1],-delr1[2]); - - // attractive term via loop over k - - for (kk = 0; kk < jnum; kk++) { - k = jlist[kk]; - k &= NEIGHMASK; - kel = type_map[type[k]]; - if (kel < 0) continue; - - 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]; - - iparam_ik = elem2param[kel][kel]; - PairParameters & q = pairParameters[iparam_ik]; - if (rsq2 > q.cutsq) continue; - r2 = sqrt(rsq2); - - fpair = (q.W)->derivative(r2); - fpair = -prefactor*fpair/r2; - - f[i][0] += delr2[0]*fpair; - f[i][1] += delr2[1]*fpair; - f[i][2] += delr2[2]*fpair; - f[k][0] -= delr2[0]*fpair; - f[k][1] -= delr2[1]*fpair; - f[k][2] -= delr2[2]*fpair; - - if (vflag_atom) v_tally2(i, k, -fpair, delr2); - } - - } else { - - // three-body interactions - // skip immediately if I-J is not within cutoff - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - jel = type_map[type[j]]; - if (jel < 0) continue; - - 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]; - - iparam_ij = elem2param[iel][jel]; - PairParameters & p = pairParameters[iparam_ij]; - if (rsq1 > p.cutsq) continue; - r1 = sqrt(rsq1); - - // accumulate bondorder zeta for each i-j interaction via loop over k - - zeta_ij = 0.0; - - for (kk = 0; kk < jnum; kk++) { - if (jj == kk) continue; - k = jlist[kk]; - k &= NEIGHMASK; - kel = type_map[type[k]]; - if (kel < 0) continue; - - 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]; - - iparam_ik = elem2param[iel][kel]; - PairParameters & q = pairParameters[iparam_ik]; - if (rsq2 > q.cutsq) continue; - r2 = sqrt(rsq2); - - costheta = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + - delr1[2]*delr2[2]) / (r1*r2); - iparam_ijk = elem3param[jel][iel][kel]; - TripletParameters & trip = tripletParameters[iparam_ijk]; - - wfac= (q.W)->value(r2); - pfac= (q.P)->value(r1-(p.xi)*r2); - gfac= (trip.G)->value(costheta); - - zeta_ij += wfac*pfac*gfac; - } - - // pairwise force due to zeta - - fa = (p.V)->value(r1); - fa_d = (p.V)->derivative(r1); - bij = (p.F)->value(zeta_ij); - bij_d = (p.F)->derivative(zeta_ij); - fpair = -0.5*bij*fa_d / r1; - prefactor = 0.5* fa * bij_d; - if (eflag) evdwl = -0.5*bij*fa; - - f[i][0] += delr1[0]*fpair; - f[i][1] += delr1[1]*fpair; - f[i][2] += delr1[2]*fpair; - f[j][0] -= delr1[0]*fpair; - f[j][1] -= delr1[1]*fpair; - f[j][2] -= delr1[2]*fpair; - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,-fpair,-delr1[0],-delr1[1],-delr1[2]); - - // attractive term via loop over k - - for (kk = 0; kk < jnum; kk++) { - if (jj == kk) continue; - k = jlist[kk]; - k &= NEIGHMASK; - kel = type_map[type[k]]; - if (kel < 0) continue; - - 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]; - - iparam_ik = elem2param[iel][kel]; - PairParameters & q = pairParameters[iparam_ik]; - if (rsq2 > q.cutsq) continue; - r2 = sqrt(rsq2); - - iparam_ijk = elem3param[jel][iel][kel]; - TripletParameters & trip = tripletParameters[iparam_ijk]; - - attractive(&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk); - - f[i][0] += fi[0]; - f[i][1] += fi[1]; - f[i][2] += fi[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 (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); - } - } - } - } - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::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"); - - type_map = new int[n+1]; -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairPolymorphic::settings(int narg, char **arg) -{ - if (narg != 0) error->all(FLERR,"Illegal pair_style command"); -} - - -/* ---------------------------------------------------------------------- - init specific to this pair style -------------------------------------------------------------------------- */ - -void PairPolymorphic::init_style() -{ - if (atom->tag_enable == 0) - error->all(FLERR,"Pair style polymorphic requires atom IDs"); - if (force->newton_pair == 0) - error->all(FLERR,"Pair style polymorphic 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 PairPolymorphic::init_one(int i, int j) -{ - if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); - - return cutmax; -} - - -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::setup() -{ - int i,j,k,n; - - memory->destroy(elem2param); - memory->create(elem2param,nelements,nelements,"pair:elem2param"); - memory->destroy(elem3param); - memory->create(elem3param,nelements,nelements,nelements,"pair:elem3param"); - - // map atom pair to parameter index, as read from potential file - n = 0; - for (i = 0; i < nelements; i++) { // note self first - elem2param[i][i] = n; - n++; - } - for (i = 0; i < nelements; i++) - for (j = i+1; j < nelements; j++) { - elem2param[i][j] = n; - elem2param[j][i] = n; - n++; - } - - // map atom triplet to parameter index - n = 0; - for (i = 0; i < nelements; i++) - for (j = 0; j < nelements; j++) - for (k = 0; k < nelements; k++) { - elem3param[i][j][k] = n; - n++; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairPolymorphic::coeff(int narg, char **arg) -{ - 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 type_map atom types to elements in potential file - int ntypes = atom->ntypes; - // type_map = atom type to element in potential file - if (type_map) { delete [] type_map; } - type_map = new int[ntypes+1]; - for (int i = 0; i < ntypes+1; i++) { - type_map[i] = -1; - } - // elements = list of requested element names (ntypes long) - char** elements = new char*[ntypes]; - for (int i = 0; i < ntypes; i++) { elements[i] = NULL; } - // parse and store - for (int i = 3; i < narg; i++) { - if (strcmp(arg[i],"NULL") != 0) { - int n = strlen(arg[i]) + 1; - elements[i-3] = new char[n]; - strcpy(elements[i-3],arg[i]); - } - } - - // read potential file and initialize potential parameters - read_file(arg[2],elements); - setup(); - - if (elements) - for (int i = 0; i < ntypes; i++) - if (elements[i]) delete [] elements[i]; - delete [] elements; - - // 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 type_mapped to elements - int count = 0; - for (int i = 1; i <= n; i++) { - for (int j = i; j <= n; j++) { - if ((type_map[i] > -1) && (type_map[j] > -1)) { - setflag[i][j] = 1; - count++; - } - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::read_file(char *file, char** elements) -{ - char line[MAXLINE],*ptr, ftype[MAXLINE]; - int n; - // open file on proc 0 - FILE *fp=NULL; - if (comm->me == 0) { - fp = force->open_potential(file); - if (fp == NULL) { - char str[128]; - sprintf(str,"Cannot open polymorphic potential file %s",file); - error->one(FLERR,str); - } - // move past comments to first data line - fgets(line,MAXLINE,fp); - while (line == strchr(line,'#')) fgets(line,MAXLINE,fp); - n = strlen(line) + 1; - } - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); - ptr = strtok(line," \t\n\r\f"); // 1st line, 1st token : nelements - nelements = atoi(ptr); // number of elements in potential file - ptr = strtok(NULL," \t\n\r\f"); // 1st line, 2nd token : indicator eta - eta = (atoi(ptr)>0) ? true:false; - if (comm->me == 0) { printf("%d elements in: %s,",nelements,file); } - - // type_map the elements in the potential file to LAMMPS atom types - for (int i = 0; i < nelements; i++) { - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token atomic number - ptr = strtok(NULL," \t\n\r\f"); // 2st token atomic mass - ptr = strtok(NULL," \t\n\r\f"); // 3st token atomic symbol - if (comm->me == 0) { printf(" %s",ptr); } - int j = 0; - for (j = 0; j < atom->ntypes; j++) { - if (elements[j] && strcmp(ptr,elements[j]) == 0) { - type_map[j+1] = i; - if (comm->me == 0) { printf("=%d ",j+1); } - break; - } - } - if (j == nelements) - error->all(FLERR,"Element not defined in potential file"); - } - if (comm->me == 0) { printf("\n"); } - - // size - npair = nelements*(nelements+1)/2; - ntriple = nelements*nelements*nelements; - pairParameters = (PairParameters*) - memory->srealloc(pairParameters,npair*sizeof(PairParameters), - "pair:pairParameters"); - tripletParameters = (TripletParameters*) - memory->srealloc(tripletParameters,ntriple*sizeof(TripletParameters), - "pair:tripletParameters"); - - // pairwise cutoffs - for (int i = 0; i < npair; i++) { - PairParameters & p = pairParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token: cutoff - p.cut = atof(ptr); - p.cutsq = p.cut*p.cut; - ptr = strtok(NULL," \t\n\r\f"); // 2nd token: indicator xi - p.xi = (atoi(ptr)>0) ? true:false; - } - - // set cutmax to max of all params - cutmax = 0.0; - for (int i = 0; i < npair; i++) { - PairParameters & p = pairParameters[i]; - if (p.cut > cutmax) cutmax = p.cut; - } - - // start reading functions - for (int i = 0; i < npair; i++) { // U - PairParameters & p = pairParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token - strcpy(ftype,ptr); - p.U = create_function(ftype,fp); - } - for (int i = 0; i < npair; i++) { // V - PairParameters & p = pairParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token - strcpy(ftype,ptr); - p.V = create_function(ftype,fp); - } - for (int i = 0; i < npair; i++) { // W - PairParameters & p = pairParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token - strcpy(ftype,ptr); - p.W = create_function(ftype,fp); - } - for (int i = 0; i < npair; i++) { // P - PairParameters & p = pairParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token - strcpy(ftype,ptr); - p.P = create_function(ftype,fp); - } - for (int i = 0; i < ntriple; i++) { // G - TripletParameters & p = tripletParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token - strcpy(ftype,ptr); - p.G = create_function(ftype,fp); - } - for (int i = 0; i < npair; i++) { // F - PairParameters & p = pairParameters[i]; - read_line(fp,line); - ptr = strtok(line," \t\n\r\f"); // 1st token - strcpy(ftype,ptr); - p.F = create_function(ftype,fp); - } - if (comm->me == 0) { fclose(fp); } -} - -/* ---------------------------------------------------------------------- */ - -C1function * PairPolymorphic::create_function(char* ftype, FILE* fp) -{ - char * ptr; - if (strcmp(ftype,"spline")==0) { // N, min, max, values - C1tabularFunction * f = new C1tabularFunction(); - ptr = strtok(NULL," \t\n\r\f"); - int n = atof(ptr); - ptr = strtok(NULL," \t\n\r\f"); - double xmin = atof(ptr); - ptr = strtok(NULL," \t\n\r\f"); - double xmax = atof(ptr); - double * table = new double[n]; - read_array(fp,n,table); - f->set_values(n,xmin,xmax,table); - delete [] table; - return f; - } - else if (strcmp(ftype,"constant") == 0) { - ptr = strtok(NULL," \t\n\r\f"); - double c = atof(ptr); - return new C1constant(c); - } - else if (strcmp(ftype,"exponential") == 0) { - ptr = strtok(NULL," \t\n\r\f"); - double c = atof(ptr); - ptr = strtok(NULL," \t\n\r\f"); - double lambda = atof(ptr); - return new C1exponential(c,lambda); - } - else if (strcmp(ftype,"sine") == 0) { - ptr = strtok(NULL," \t\n\r\f"); - double c = atof(ptr); - ptr = strtok(NULL," \t\n\r\f"); - double w = atof(ptr); - return new C1sine(c,w); - } - else if (strcmp(ftype,"cosine") == 0) { - ptr = strtok(NULL," \t\n\r\f"); - double c = atof(ptr); - ptr = strtok(NULL," \t\n\r\f"); - double w = atof(ptr); - return new C1cosine(c,w); - } - else { error->all(FLERR,"unknown function type"); } - return NULL; -} -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::read_line(FILE *fp, char *line) -{ - int n = 0; - if (comm->me == 0) { - fgets(line,MAXLINE,fp); - n = strlen(line) + 1; - } - MPI_Bcast(&n,1,MPI_INT,0,world); - MPI_Bcast(line,n,MPI_CHAR,0,world); -} -void PairPolymorphic::read_array(FILE *fp, int n, double *list) -{ - if (comm->me == 0) { - char *ptr; - char line[MAXLINE]; - int i = 0; - while (i < n) { - fgets(line,MAXLINE,fp); - ptr = strtok(line," \t\n\r\f"); - list[i++] = atof(ptr); - while ((ptr = strtok(NULL," \t\n\r\f"))) list[i++] = atof(ptr); - } - } - MPI_Bcast(list,n,MPI_DOUBLE,0,world); -} - - -/* ---------------------------------------------------------------------- - attractive term -------------------------------------------------------------------------- */ - -void PairPolymorphic::attractive(PairParameters *p, TripletParameters *trip, - double prefactor, double rij, double rik, - double *delrij, double *delrik, - double *fi, double *fj, double *fk) -{ - double rij_hat[3],rik_hat[3]; - double rijinv,rikinv; - - rijinv = 1.0/rij; - vec3_scale(rijinv,delrij,rij_hat); - - rikinv = 1.0/rik; - vec3_scale(rikinv,delrik,rik_hat); - - ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,p,trip); -} - -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::ters_zetaterm_d(double prefactor, - double *rij_hat, double rij, - double *rik_hat, double rik, - double *dri, double *drj, double *drk, - PairParameters *p, TripletParameters *trip) -{ - double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta; - double dcosdri[3],dcosdrj[3],dcosdrk[3]; - - cos_theta = vec3_dot(rij_hat,rik_hat); - - fc = (p->W)->value(rik); - dfc = (p->W)->derivative(rik); - ex_delr = (p->P)->value(rij-(p->xi)*rik); - ex_delr_d = (p->P)->derivative(rij-(p->xi)*rik); - gijk = (trip->G)->value(cos_theta); - gijk_d = (trip->G)->derivative(cos_theta); - - costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); - - // compute the derivative wrt Ri - // dri = -dfc*gijk*ex_delr*rik_hat; - // dri += fc*gijk_d*ex_delr*dcosdri; - // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat); - - vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri); - vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri); - vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri); - vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri); - vec3_scale(prefactor,dri,dri); - - // compute the derivative wrt Rj - // drj = fc*gijk_d*ex_delr*dcosdrj; - // drj += fc*gijk*ex_delr_d*rij_hat; - - vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj); - vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj); - vec3_scale(prefactor,drj,drj); - - // compute the derivative wrt Rk - // drk = dfc*gijk*ex_delr*rik_hat; - // drk += fc*gijk_d*ex_delr*dcosdrk; - // drk += -fc*gijk*ex_delr_d*rik_hat; - - vec3_scale(dfc*gijk*ex_delr,rik_hat,drk); - vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk); - vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk); - vec3_scale(prefactor,drk,drk); -} - -/* ---------------------------------------------------------------------- */ - -void PairPolymorphic::costheta_d(double *rij_hat, double rij, - double *rik_hat, double rik, - double *dri, double *drj, double *drk) -{ - // first element is devative wrt Ri, second wrt Rj, third wrt Rk - - double cos_theta = vec3_dot(rij_hat,rik_hat); - - vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj); - vec3_scale(1.0/rij,drj,drj); - vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk); - vec3_scale(1.0/rik,drk,drk); - vec3_add(drj,drk,dri); - vec3_scale(-1.0,dri,dri); -} - +/* ---------------------------------------------------------------------- + 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: Reese Jones, Xiaowang Zhou (SNL) + This modifies from pair_tersoff.cpp by Aidan Thompson (SNL) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_polymorphic.h" +#include "atom.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "error.h" + +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define MAXLINE 1024 +#define DELTA 4 + +/* ---------------------------------------------------------------------- */ + +PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + one_coeff = 1; + + nelements = 0; + elements = NULL; + pairParameters = NULL; + tripletParameters = NULL; + elem2param = NULL; + elem3param = NULL; + epsilon = 0.0; + neighsize = 0; + firstneighV = NULL; + firstneighW = NULL; + firstneighW1 = NULL; + delxV = NULL; + delyV = NULL; + delzV = NULL; + drV = NULL; + delxW = NULL; + delyW = NULL; + delzW = NULL; + drW = NULL; +} + +/* ---------------------------------------------------------------------- + check if allocated, since class can be destructed when incomplete +------------------------------------------------------------------------- */ + +PairPolymorphic::~PairPolymorphic() +{ + if (elements) + for (int i = 0; i < nelements; i++) delete [] elements[i]; + delete [] elements; + delete [] match; + memory->destroy(pairParameters); + memory->destroy(tripletParameters); + memory->destroy(elem2param); + memory->destroy(elem3param); + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + delete [] map; + delete [] firstneighV; + delete [] firstneighW; + delete [] firstneighW1; + delete [] delxV; + delete [] delyV; + delete [] delzV; + delete [] drV; + delete [] delxW; + delete [] delyW; + delete [] delzW; + delete [] drW; + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::compute(int eflag, int vflag) +{ + int i,j,k,ii,jj,kk,kk1,inum,jnum; + int itag,jtag,itype,jtype,ktype; + int iparam_ii,iparam_jj,iparam_kk,iparam_ij,iparam_ik,iparam_ijk; + double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; + double rsq,rsq1,rsq2,r0,r1,r2; + double delr1[3],delr2[3],fi[3],fj[3],fk[3]; + double zeta_ij,prefactor,wfac,pfac,gfac,fa,fa_d,bij,bij_d; + double costheta; + int *ilist,*jlist,*numneigh,**firstneigh; + double emb; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = vflag_atom = 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]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + if (neighsize < jnum) { + delete [] firstneighV; + delete [] delxV; + delete [] delyV; + delete [] delzV; + delete [] drV; + delete [] firstneighW; + delete [] delxW; + delete [] delyW; + delete [] delzW; + delete [] drW; + delete [] firstneighW1; + neighsize = jnum + 20; + firstneighV = new int[neighsize]; + delxV = new double[neighsize]; + delyV = new double[neighsize]; + delzV = new double[neighsize]; + drV = new double[neighsize]; + firstneighW = new int[neighsize]; + delxW = new double[neighsize]; + delyW = new double[neighsize]; + delzW = new double[neighsize]; + drW = new double[neighsize]; + firstneighW1 = new int[neighsize]; + } + + if (eta) { + iparam_ii = elem2param[itype][itype]; + PairParameters & p = pairParameters[iparam_ii]; + emb = (p.F)->get_vmax(); + } + + numneighV = -1; + numneighW = -1; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + 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; + if (rsq >= cutmaxsq) continue; + r0 = sqrt(rsq); + + iparam_ij = elem2param[itype][jtype]; + PairParameters & p = pairParameters[iparam_ij]; + +// do not include the neighbor if get_vmax() <= epsilon because the function is near zero + if (eta) { + if (emb > epsilon) { + iparam_jj = elem2param[jtype][jtype]; + PairParameters & q = pairParameters[iparam_jj]; + if (rsq < (q.W)->get_xmaxsq() && (q.W)->get_vmax() > epsilon) { + numneighW = numneighW + 1; + firstneighW[numneighW] = j; + delxW[numneighW] = delx; + delyW[numneighW] = dely; + delzW[numneighW] = delz; + drW[numneighW] = r0; + } + } + } else { + if ((p.F)->get_vmax() > epsilon) { + if (rsq < (p.V)->get_xmaxsq() && (p.V)->get_vmax() > epsilon) { + numneighV = numneighV + 1; + firstneighV[numneighV] = j; + delxV[numneighV] = delx; + delyV[numneighV] = dely; + delzV[numneighV] = delz; + drV[numneighV] = r0; + } + if (rsq < (p.W)->get_xmaxsq() && (p.W)->get_vmax() > epsilon) { + numneighW = numneighW + 1; + firstneighW[numneighW] = j; + delxW[numneighW] = delx; + delyW[numneighW] = dely; + delzW[numneighW] = delz; + drW[numneighW] = r0; + } + } + } + + // two-body interactions, skip half of them + + 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] < x[i][2]) 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; + } + + if (rsq >= (p.U)->get_xmaxsq() || (p.U)->get_vmax() <= epsilon) continue; + (p.U)->value(r0,evdwl,eflag,fpair,1); + fpair = -fpair/r0; + + 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); + } + + if (eta) { + + if (emb > epsilon) { + + iparam_ii = elem2param[itype][itype]; + PairParameters & p = pairParameters[iparam_ii]; + + // accumulate bondorder zeta for each i-j interaction via loop over k + + zeta_ij = 0.0; + + for (kk = 0; kk <= numneighW; kk++) { + k = firstneighW[kk]; + ktype = map[type[k]]; + + iparam_kk = elem2param[ktype][ktype]; + PairParameters & q = pairParameters[iparam_kk]; + + (q.W)->value(drW[kk],wfac,1,fpair,0); + + zeta_ij += wfac; + } + + // pairwise force due to zeta + + (p.F)->value(zeta_ij,bij,1,bij_d,1); + + prefactor = 0.5* bij_d; + if (eflag) evdwl = -0.5*bij; + + if (evflag) ev_tally(i,i,nlocal,newton_pair,evdwl,0.0,0.0,delx,dely,delz); + + // attractive term via loop over k + + for (kk = 0; kk <= numneighW; kk++) { + k = firstneighW[kk]; + ktype = map[type[k]]; + + delr2[0] = -delxW[kk]; + delr2[1] = -delyW[kk]; + delr2[2] = -delzW[kk]; + + iparam_kk = elem2param[ktype][ktype]; + PairParameters & q = pairParameters[iparam_kk]; + + (q.W)->value(drW[kk],wfac,0,fpair,1); + fpair = -prefactor*fpair/drW[kk]; + + f[i][0] += delr2[0]*fpair; + f[i][1] += delr2[1]*fpair; + f[i][2] += delr2[2]*fpair; + f[k][0] -= delr2[0]*fpair; + f[k][1] -= delr2[1]*fpair; + f[k][2] -= delr2[2]*fpair; + + if (vflag_atom) v_tally2(i, k, -fpair, delr2); + } + } + + } else { + + for (jj = 0; jj <= numneighV; jj++) { + j = firstneighV[jj]; + jtype = map[type[j]]; + + iparam_ij = elem2param[itype][jtype]; + PairParameters & p = pairParameters[iparam_ij]; + + delr1[0] = -delxV[jj]; + delr1[1] = -delyV[jj]; + delr1[2] = -delzV[jj]; + r1 = drV[jj]; + + // accumulate bondorder zeta for each i-j interaction via loop over k + + zeta_ij = 0.0; + + numneighW1 = -1; + for (kk = 0; kk <= numneighW; kk++) { + k = firstneighW[kk]; + if (j == k) continue; + ktype = map[type[k]]; + iparam_ijk = elem3param[jtype][itype][ktype]; + TripletParameters & trip = tripletParameters[iparam_ijk]; + if ((trip.G)->get_vmax() <= epsilon) continue; + + numneighW1 = numneighW1 + 1; + firstneighW1[numneighW1] = kk; + + delr2[0] = -delxW[kk]; + delr2[1] = -delyW[kk]; + delr2[2] = -delzW[kk]; + r2 = drW[kk]; + + costheta = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + + delr1[2]*delr2[2]) / (r1*r2); + + iparam_ik = elem2param[itype][ktype]; + PairParameters & q = pairParameters[iparam_ik]; + + (q.W)->value(r2,wfac,1,fpair,0); + (q.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0); + (trip.G)->value(costheta,gfac,1,fpair,0); + + zeta_ij += wfac*pfac*gfac; + } + + // pairwise force due to zeta + + (p.V)->value(r1,fa,1,fa_d,1); + (p.F)->value(zeta_ij,bij,1,bij_d,1); + fpair = -0.5*bij*fa_d / r1; + prefactor = 0.5* fa * bij_d; + if (eflag) evdwl = -0.5*bij*fa; + + f[i][0] += delr1[0]*fpair; + f[i][1] += delr1[1]*fpair; + f[i][2] += delr1[2]*fpair; + f[j][0] -= delr1[0]*fpair; + f[j][1] -= delr1[1]*fpair; + f[j][2] -= delr1[2]*fpair; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,-fpair,-delr1[0],-delr1[1],-delr1[2]); + + // attractive term via loop over k + + for (kk1 = 0; kk1 <= numneighW1; kk1++) { + kk = firstneighW1[kk1]; + k = firstneighW[kk]; + ktype = map[type[k]]; + iparam_ijk = elem3param[jtype][itype][ktype]; + TripletParameters & trip = tripletParameters[iparam_ijk]; + + delr2[0] = -delxW[kk]; + delr2[1] = -delyW[kk]; + delr2[2] = -delzW[kk]; + r2 = drW[kk]; + + iparam_ik = elem2param[itype][ktype]; + PairParameters & q = pairParameters[iparam_ik]; + + attractive(&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk); + + f[i][0] += fi[0]; + f[i][1] += fi[1]; + f[i][2] += fi[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 (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); + } + } + } + } + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::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]; + + neighsize = 40; + firstneighV = new int[neighsize]; + delxV = new double[neighsize]; + delyV = new double[neighsize]; + delzV = new double[neighsize]; + drV = new double[neighsize]; + firstneighW = new int[neighsize]; + delxW = new double[neighsize]; + delyW = new double[neighsize]; + delzW = new double[neighsize]; + drW = new double[neighsize]; + firstneighW1 = new int[neighsize]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairPolymorphic::settings(int narg, char **arg) +{ + if (narg != 0) error->all(FLERR,"Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairPolymorphic::coeff(int narg, char **arg) +{ + int i,j,n; + + if (!allocated) allocate(); + + if (narg == 4 + atom->ntypes) { + narg--; + epsilon = atof(arg[narg]); + } else 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 PairPolymorphic::init_style() +{ + if (atom->tag_enable == 0) + error->all(FLERR,"Pair style eqfree requires atom IDs"); + if (force->newton_pair == 0) + error->all(FLERR,"Pair style eqfree 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 PairPolymorphic::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + return cutmax; +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::read_file(char *file) +{ + char line[MAXLINE],*ptr; + int n; + + // open file on proc 0 + FILE *fp=NULL; + if (comm->me == 0) { + fp = fopen(file,"r"); + if (fp == NULL) { + char str[128]; + sprintf(str,"Cannot open eqfree potential file %s",file); + error->one(FLERR,str); + } + // move past comments to first data line + fgets(line,MAXLINE,fp); + while (line == strchr(line,'#')) fgets(line,MAXLINE,fp); + n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + ptr = strtok(line," \t\n\r\f"); // 1st line, 1st token + int ntypes = atoi(ptr); + if (ntypes != nelements) + error->all(FLERR,"Incorrect number of elements in potential file"); + match = new int[nelements]; + ptr = strtok(NULL," \t\n\r\f"); // 1st line, 2nd token + eta = (atoi(ptr)>0) ? true:false; + + // map the elements in the potential file to LAMMPS atom types + for (int i = 0; i < nelements; i++) { + if (comm->me == 0) { + fgets(line,MAXLINE,fp); + n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + ptr = strtok(line," \t\n\r\f"); // 1st token + ptr = strtok(NULL," \t\n\r\f"); // 2st token + ptr = strtok(NULL," \t\n\r\f"); // 3st token + int j; + for (j = 0; j < nelements; j++) { + if (strcmp(ptr,elements[j]) == 0) break; + } + if (j == nelements) + error->all(FLERR,"Element not defined in potential file"); + match[i] = j; + } + // sizes + if (comm->me == 0) { + fgets(line,MAXLINE,fp); + n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + ptr = strtok(line," \t\n\r\f"); // 1st token + nr = atoi(ptr); + ptr = strtok(NULL," \t\n\r\f"); // 2nd token + ng = atoi(ptr); + ptr = strtok(NULL," \t\n\r\f"); // 3rd token + nx = atoi(ptr); + ptr = strtok(NULL," \t\n\r\f"); // 4th token + maxX = atof(ptr); + + npair = nelements*(nelements+1)/2; + ntriple = nelements*nelements*nelements; + pairParameters = (PairParameters*) + memory->srealloc(pairParameters,npair*sizeof(PairParameters), + "pair:pairParameters"); + tripletParameters = (TripletParameters*) + memory->srealloc(tripletParameters,ntriple*sizeof(TripletParameters), + "pair:tripletParameters"); + + // cutoffs + for (int i = 0; i < npair; i++) { + PairParameters & p = pairParameters[i]; + if (comm->me == 0) { + fgets(line,MAXLINE,fp); + n = strlen(line) + 1; + } + MPI_Bcast(&n,1,MPI_INT,0,world); + MPI_Bcast(line,n,MPI_CHAR,0,world); + ptr = strtok(line," \t\n\r\f"); // 1st token + p.cut = atof(ptr); + p.cutsq = p.cut*p.cut; + ptr = strtok(NULL," \t\n\r\f"); // 2nd token + p.xi = (atoi(ptr)>0) ? true:false; + } + + // set cutmax to max of all params + cutmax = 0.0; + for (int i = 0; i < npair; i++) { + PairParameters & p = pairParameters[i]; + if (p.cut > cutmax) cutmax = p.cut; + } + cutmaxsq = cutmax*cutmax; + + // start reading tabular functions + double * singletable = new double[nr]; + for (int i = 0; i < npair; i++) { // U + PairParameters & p = pairParameters[i]; + if (comm->me == 0) { + grab(fp,nr,singletable); + } + MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world); + p.U = new tabularFunction(nr,0.0,p.cut); + (p.U)->set_values(nr,0.0,p.cut,singletable,epsilon); + } + for (int i = 0; i < npair; i++) { // V + PairParameters & p = pairParameters[i]; + if (comm->me == 0) { + grab(fp,nr,singletable); + } + MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world); + p.V = new tabularFunction(nr,0.0,p.cut); + (p.V)->set_values(nr,0.0,p.cut,singletable,epsilon); + } + for (int i = 0; i < npair; i++) { // W + PairParameters & p = pairParameters[i]; + if (comm->me == 0) { + grab(fp,nr,singletable); + } + MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world); + p.W = new tabularFunction(nr,0.0,p.cut); + (p.W)->set_values(nr,0.0,p.cut,singletable,epsilon); + } + for (int i = 0; i < npair; i++) { // P + PairParameters & p = pairParameters[i]; + if (comm->me == 0) { + grab(fp,nr,singletable); + } + MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world); + p.P = new tabularFunction(nr,-cutmax,cutmax); + (p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon); + } + delete singletable; + singletable = new double[ng]; + for (int i = 0; i < ntriple; i++) { // G + TripletParameters & p = tripletParameters[i]; + if (comm->me == 0) { + grab(fp,ng,singletable); + } + MPI_Bcast(singletable,ng,MPI_DOUBLE,0,world); + p.G = new tabularFunction(ng,-1.0,1.0); + (p.G)->set_values(ng,-1.0,1.0,singletable,epsilon); + } + delete singletable; + singletable = new double[nx]; + for (int i = 0; i < npair; i++) { // F + PairParameters & p = pairParameters[i]; + if (comm->me == 0) { + grab(fp,nx,singletable); + } + MPI_Bcast(singletable,nx,MPI_DOUBLE,0,world); + p.F = new tabularFunction(nx,0.0,maxX); + (p.F)->set_values(nx,0.0,maxX,singletable,epsilon); + } + delete singletable; + if (comm->me == 0) { + fclose(fp); + } + +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::setup() +{ + int i,j,k,n; + + memory->destroy(elem2param); + memory->create(elem2param,nelements,nelements,"pair:elem2param"); + memory->destroy(elem3param); + memory->create(elem3param,nelements,nelements,nelements,"pair:elem3param"); + + // map atom pair to parameter index + + n = 0; + for (i = 0; i < nelements; i++) { + elem2param[match[i]][match[i]] = n; + n++; + } + for (i = 0; i < nelements-1; i++) { + for (j = i+1; j < nelements; j++) { + elem2param[match[i]][match[j]] = n; + elem2param[match[j]][match[i]] = n; + n++; + } + } + + // map atom triplet to parameter index + + n = 0; + for (i = 0; i < nelements; i++) + for (j = 0; j < nelements; j++) + for (k = 0; k < nelements; k++) { + elem3param[match[i]][match[j]][match[k]] = n; + n++; + } + +// for debugging, call write_tables() to check the tabular functions +// if (comm->me == 0) { +// write_tables(51); +// error->all(FLERR,"Test potential tables"); +// } +} + +/* ---------------------------------------------------------------------- + attractive term +------------------------------------------------------------------------- */ + +void PairPolymorphic::attractive(PairParameters *p, TripletParameters *trip, + double prefactor, double rij, double rik, + double *delrij, double *delrik, + double *fi, double *fj, double *fk) +{ + double rij_hat[3],rik_hat[3]; + double rijinv,rikinv; + + rijinv = 1.0/rij; + vec3_scale(rijinv,delrij,rij_hat); + + rikinv = 1.0/rik; + vec3_scale(rikinv,delrik,rik_hat); + + ters_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,p,trip); +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::ters_zetaterm_d(double prefactor, + double *rij_hat, double rij, + double *rik_hat, double rik, + double *dri, double *drj, double *drk, + PairParameters *p, TripletParameters *trip) +{ + double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta; + double dcosdri[3],dcosdrj[3],dcosdrk[3]; + + cos_theta = vec3_dot(rij_hat,rik_hat); + + (p->W)->value(rik,fc,1,dfc,1); + (p->P)->value(rij-(p->xi)*rik,ex_delr,1,ex_delr_d,1); + (trip->G)->value(cos_theta,gijk,1,gijk_d,1); + + costheta_d(rij_hat,rij,rik_hat,rik,dcosdri,dcosdrj,dcosdrk); + + // compute the derivative wrt Ri + // dri = -dfc*gijk*ex_delr*rik_hat; + // dri += fc*gijk_d*ex_delr*dcosdri; + // dri += fc*gijk*ex_delr_d*(rik_hat - rij_hat); + + vec3_scale(-dfc*gijk*ex_delr,rik_hat,dri); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdri,dri,dri); + vec3_scaleadd(fc*gijk*ex_delr_d,rik_hat,dri,dri); + vec3_scaleadd(-fc*gijk*ex_delr_d,rij_hat,dri,dri); + vec3_scale(prefactor,dri,dri); + + // compute the derivative wrt Rj + // drj = fc*gijk_d*ex_delr*dcosdrj; + // drj += fc*gijk*ex_delr_d*rij_hat; + + vec3_scale(fc*gijk_d*ex_delr,dcosdrj,drj); + vec3_scaleadd(fc*gijk*ex_delr_d,rij_hat,drj,drj); + vec3_scale(prefactor,drj,drj); + + // compute the derivative wrt Rk + // drk = dfc*gijk*ex_delr*rik_hat; + // drk += fc*gijk_d*ex_delr*dcosdrk; + // drk += -fc*gijk*ex_delr_d*rik_hat; + + vec3_scale(dfc*gijk*ex_delr,rik_hat,drk); + vec3_scaleadd(fc*gijk_d*ex_delr,dcosdrk,drk,drk); + vec3_scaleadd(-fc*gijk*ex_delr_d,rik_hat,drk,drk); + vec3_scale(prefactor,drk,drk); +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::costheta_d(double *rij_hat, double rij, + double *rik_hat, double rik, + double *dri, double *drj, double *drk) +{ + // first element is devative wrt Ri, second wrt Rj, third wrt Rk + + double cos_theta = vec3_dot(rij_hat,rik_hat); + + vec3_scaleadd(-cos_theta,rij_hat,rik_hat,drj); + vec3_scale(1.0/rij,drj,drj); + vec3_scaleadd(-cos_theta,rik_hat,rij_hat,drk); + vec3_scale(1.0/rik,drk,drk); + vec3_add(drj,drk,dri); + vec3_scale(-1.0,dri,dri); +} + +/* ---------------------------------------------------------------------- + * grab n values from file fp and put them in list + * values can be several to a line + * only called by proc 0 + * ------------------------------------------------------------------------- */ + +void PairPolymorphic::grab(FILE *fp, int n, double *list) +{ + char *ptr; + char line[MAXLINE]; + + int i = 0; + while (i < n) { + fgets(line,MAXLINE,fp); + ptr = strtok(line," \t\n\r\f"); + list[i++] = atof(ptr); + while (ptr = strtok(NULL," \t\n\r\f")) list[i++] = atof(ptr); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPolymorphic::write_tables(int npts) +{ + char tag[6] = ""; + if (comm->me != 0) sprintf(tag,"%d",comm->me); + FILE* fp = NULL; + double xmin,xmax,x,uf,vf,wf,pf,gf,ff,ufp,vfp,wfp,pfp,gfp,ffp; + char line[MAXLINE]; + for (int i = 0; i < nelements; i++) { + for (int j = 0; j < nelements; j++) { + strcpy(line,elements[i]); + strcat(line,elements[j]); + strcat(line,"_UVW"); + strcat(line,tag); + fp = fopen(line, "w"); + int iparam_ij = elem2param[i][j]; + PairParameters & pair = pairParameters[iparam_ij]; + xmin = (pair.U)->get_xmin(); + xmax = (pair.U)->get_xmax(); + double xl = xmax - xmin; + xmin = xmin - 0.5*xl; + xmax = xmax + 0.5*xl; + for (int k = 0; k < npts; k++) { + x = xmin + (xmax-xmin) * k / (npts-1); + (pair.U)->value(x, uf, 1, ufp, 1); + (pair.V)->value(x, vf, 1, vfp, 1); + (pair.W)->value(x, wf, 1, wfp, 1); + fprintf(fp,"%12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f \n",x,uf,vf,wf,ufp,vfp,wfp); + } + fclose(fp); + } + } + for (int i = 0; i < nelements; i++) { + for (int j = 0; j < nelements; j++) { + strcpy(line,elements[i]); + strcat(line,elements[j]); + strcat(line,"_P"); + strcat(line,tag); + fp = fopen(line, "w"); + int iparam_ij = elem2param[i][j]; + PairParameters & pair = pairParameters[iparam_ij]; + xmin = (pair.P)->get_xmin(); + xmax = (pair.P)->get_xmax(); + double xl = xmax - xmin; + xmin = xmin - 0.5*xl; + xmax = xmax + 0.5*xl; + for (int k = 0; k < npts; k++) { + x = xmin + (xmax-xmin) * k / (npts-1); + (pair.P)->value(x, pf, 1, pfp, 1); + fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp); + } + fclose(fp); + } + } + for (int i = 0; i < nelements; i++) { + for (int j = 0; j < nelements; j++) { + for (int k = 0; k < nelements; k++) { + strcpy(line,elements[i]); + strcat(line,elements[j]); + strcat(line,elements[k]); + strcat(line,"_G"); + strcat(line,tag); + fp = fopen(line, "w"); + int iparam_ij = elem3param[i][j][k]; + TripletParameters & pair = tripletParameters[iparam_ij]; + xmin = (pair.G)->get_xmin(); + xmax = (pair.G)->get_xmax(); + for (int n = 0; n < npts; n++) { + x = xmin + (xmax-xmin) * n / (npts-1); + (pair.G)->value(x, gf, 1, gfp, 1); + fprintf(fp,"%12.4f %12.4f %12.4f \n",x,gf,gfp); + } + fclose(fp); + } + } + } + for (int i = 0; i < nelements; i++) { + for (int j = 0; j < nelements; j++) { + strcpy(line,elements[i]); + strcat(line,elements[j]); + strcat(line,"_F"); + strcat(line,tag); + fp = fopen(line, "w"); + int iparam_ij = elem2param[i][j]; + PairParameters & pair = pairParameters[iparam_ij]; + xmin = (pair.F)->get_xmin(); + xmax = (pair.F)->get_xmax(); + double xl = xmax - xmin; + xmin = xmin - 0.5*xl; + xmax = xmax + 0.5*xl; + for (int k = 0; k < npts; k++) { + x = xmin + (xmax-xmin) * k / (npts-1); + (pair.F)->value(x, ff, 1, ffp, 1); + fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp); + } + fclose(fp); + } + } + +} + diff --git a/src/MANYBODY/pair_polymorphic.h b/src/MANYBODY/pair_polymorphic.h index 901fdaecc..6db8d9a07 100755 --- a/src/MANYBODY/pair_polymorphic.h +++ b/src/MANYBODY/pair_polymorphic.h @@ -1,315 +1,319 @@ -/* ---------------------------------------------------------------------- - 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(polymorphic,PairPolymorphic) - -#else - -#ifndef LMP_PAIR_POLYMORPHIC_H -#define LMP_PAIR_POLYMORPHIC_H - -#include "pair.h" - -namespace LAMMPS_NS { - -//=========================================== -class C1function { - public: - C1function(){}; - virtual ~C1function() {}; - virtual double value(double x)=0; - virtual double derivative(double x)=0; -}; -//=========================================== -class C1constant : public C1function { - public: - C1constant(double C): C_(C) {}; - C1constant(): C_(0) {}; - virtual double value(double x) { return C_; } - virtual double derivative(double x) { return 0.; } - protected: - double C_; -}; -//=========================================== -class C1exponential : public C1function { - public: - C1exponential(double C, double lambda): C_(C),lambda_(lambda) {}; - C1exponential(): C_(0),lambda_(0) {}; - virtual double value(double x) { return C_*exp(lambda_*x); } - virtual double derivative(double x) { return lambda_*C_*exp(lambda_*x); } - protected: - double C_,lambda_; -}; -//=========================================== -class C1sine : public C1function { - public: - C1sine(double C, double lambda): C_(C),lambda_(lambda) {}; - C1sine(): C_(0),lambda_(0) {}; - virtual double value(double x) { return C_*sin(lambda_*x); } - virtual double derivative(double x) { return lambda_*C_*cos(lambda_*x); } - protected: - double C_,lambda_; -}; -//=========================================== -class C1cosine : public C1function { - public: - C1cosine(double C, double lambda): C_(C),lambda_(lambda) {}; - C1cosine(): C_(0),lambda_(0) {}; - virtual double value(double x) { return C_*cos(lambda_*x); } - virtual double derivative(double x) { return -lambda_*C_*sin(lambda_*x); } - protected: - double C_,lambda_; -}; -//=========================================== -class C1tabularFunction : public C1function { - public: - C1tabularFunction():size(0),xmin(0),xmax(0) { - resize(0); - } - C1tabularFunction(int n):size(0),xmin(0),xmax(0) { - resize(n); - } - C1tabularFunction(int n, double x1, double x2):size(0),xmin(x1),xmax(x2) { - resize(n); - } - virtual ~C1tabularFunction() { - resize(0); - } - virtual double value(double x) { - double y,dy; - this->value(x,y,1,dy,0); - return y; - } - virtual double derivative(double x) { - double y,dy; - this->value(x,y,0,dy,1); - return dy; - } - void value(double x, double &y, int ny, double &y1, int ny1) - { - double ps = (x - xmin) * rdx; - int ks = ps; - if (ks > size-2) ks = size-2; - ps = ps - ks; - if (ps > 1.0) ps = 1.0; - if (ny) y = ((ys3[ks]*ps + ys2[ks])*ps + ys1[ks])*ps + ys[ks]; - if (ny1) y1 = (ys6[ks]*ps + ys5[ks])*ps + ys4[ks]; - } - void set_xrange(double x1, double x2) { - xmin = x1; - xmax = x2; - } - void set_values(int n, double x1, double x2, double * values) - { - resize(n); - xmin = x1; - xmax = x2; - memcpy(ys,values,n*sizeof(double)); - initialize(); - } - - void print_value() - { - printf("%d %f %f %f \n",size,xmin,xmax,rdx); - printf(" \n"); - for (int i = 0; i < size; i++) { - printf("%f %f \n",xs[i],ys[i]); - } - } - - double get_xmin() { - return xmin; - } - double get_xmax() { - return xmax; - } - double get_rdx() { - return rdx; - } - - protected: - void resize(int n) { - if (size == 0) { - xs = NULL; - ys = NULL; - ys1 = NULL; - ys2 = NULL; - ys3 = NULL; - ys4 = NULL; - ys5 = NULL; - ys6 = NULL; - } - if (n == 0) { - if (xs) delete [] xs; - if (ys) delete [] ys; - if (ys1) delete [] ys1; - if (ys2) delete [] ys2; - if (ys3) delete [] ys3; - if (ys4) delete [] ys4; - if (ys5) delete [] ys5; - if (ys6) delete [] ys6; - size = 0; - } - else if (n != size) { - size = n; - if (xs) delete [] xs; - xs = new double[n]; - if (ys) delete [] ys; - ys = new double[n]; - if (ys1) delete [] ys1; - ys1 = new double[n]; - if (ys2) delete [] ys2; - ys2 = new double[n]; - if (ys3) delete [] ys3; - ys3 = new double[n]; - if (ys4) delete [] ys4; - ys4 = new double[n]; - if (ys5) delete [] ys5; - ys5 = new double[n]; - if (ys6) delete [] ys6; - ys6 = new double[n]; - } - } - void initialize() { - int n = size; - rdx = (xmax-xmin)/(n-1.0); - for (int i = 0; i < n; i++) { - xs[i] = xmin+i*rdx; - } - rdx = 1.0 / rdx; - ys1[0] = ys[1] - ys[0]; - ys1[1] = 0.5 * (ys[2] - ys[0]); - ys1[n-2] = 0.5 * (ys[n-1] - ys[n-3]); - ys1[n-1] = ys[n-1] - ys[n-2]; - for (int i = 2; i < n-2; i++) { - ys1[i]=((ys[i-2]-ys[i+2])+ 8.0*(ys[i+1]-ys[i-1]))/12.0; - } - for (int i = 0; i < n-1; i++) { - ys2[i]=3.0*(ys[i+1]-ys[i])-2.0*ys1[i]-ys1[i+1]; - ys3[i]=ys1[i]+ys1[i+1]-2.0*(ys[i+1]-ys[i]); - } - ys2[n-1]=0.0; - ys3[n-1]=0.0; - for (int i = 0; i < n; i++) { - ys4[i]=ys1[i]*rdx; - ys5[i]=2.0*ys2[i]*rdx; - ys6[i]=3.0*ys3[i]*rdx; - } - } - int size; - double xmin,xmax,rdx; - double * ys, * ys1, * ys2, * ys3, * ys4, * ys5, * ys6; - double * xs; -}; - -//=========================================== -class PairPolymorphic : public Pair { - public: - PairPolymorphic(class LAMMPS *); - virtual ~PairPolymorphic(); - virtual void compute(int, int); - void settings(int, char **); - void coeff(int, char **); - void init_style(); - double init_one(int, int); - - protected: - struct PairParameters { - double cut; - double cutsq; - bool xi; // "indicator" - class C1function * U; - class C1function * V; - class C1function * W; - class C1function * P; - class C1function * F; - PairParameters() { - cut = 0.0; - cutsq = 0.0; - xi = true; - U = NULL; - V = NULL; - W = NULL; - P = NULL; - F = NULL; - }; - }; - struct TripletParameters { - class C1function * G; - TripletParameters() { - G = NULL; - }; - }; - - bool eta; // global indicator - int nx,nr,ng; // table sizes - double maxX; - - // parameter sets - PairParameters * pairParameters; // for I-J interaction - TripletParameters * tripletParameters; // for I-J-K interaction - - char **elements; // names of unique elements - int **elem2param; // map: element pairs to parameters - int ***elem3param; // map: element triplets to parameters - int *type_map; // mapping from atom types to elements - double cutmax; // max cutoff for all elements - int nelements; // # of unique elements - int npair,ntriple; - - void allocate(); - void read_line(FILE *, char *); - void read_array(FILE *, int, double *); - virtual void read_file(char *, char**); - class C1function * create_function(char *, FILE *); - - void setup(); - - void attractive(PairParameters *, TripletParameters *, double, double, - double, double *, double *, double *, double *, double *); - - void ters_zetaterm_d(double, double *, double, double *, double, double *, - double *, double *, PairParameters *, TripletParameters *); - void costheta_d(double *, double, double *, double, - double *, double *, double *); - - // inlined functions for efficiency - - inline double vec3_dot(const double x[3], const double y[3]) const { - return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; - } - - inline void vec3_add(const double x[3], const double y[3], - double * const z) const { - z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2]; - } - - inline void vec3_scale(const double k, const double x[3], - double y[3]) const { - y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2]; - } - - inline void vec3_scaleadd(const double k, const double x[3], - const double y[3], double * const z) const { - z[0] = k*x[0]+y[0]; - z[1] = k*x[1]+y[1]; - z[2] = k*x[2]+y[2]; - } -}; - -} - -#endif -#endif +/* ---------------------------------------------------------------------- + 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(polymorphic,PairPolymorphic) + +#else + +#ifndef LMP_PAIR_POLYMORPHIC_H +#define LMP_PAIR_POLYMORPHIC_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairPolymorphic : public Pair { + + public: + + PairPolymorphic(class LAMMPS *); + virtual ~PairPolymorphic(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + void init_style(); + double init_one(int, int); + + protected: + + class tabularFunction { + + public: + + tabularFunction() { + size = 0; + xmin = 0.0; + xmax = 0.0; + xmaxsq = xmax*xmax; + vmax = 0.0; + xs = NULL; + ys = NULL; + ys1 = NULL; + ys2 = NULL; + ys3 = NULL; + ys4 = NULL; + ys5 = NULL; + ys6 = NULL; + } + tabularFunction(int n) { + size = n; + xmin = 0.0; + xmax = 0.0; + xmaxsq = xmax*xmax; + xs = new double[n]; + ys = new double[n]; + ys1 = new double[n]; + ys2 = new double[n]; + ys3 = new double[n]; + ys4 = new double[n]; + ys5 = new double[n]; + ys6 = new double[n]; + } + tabularFunction(int n, double x1, double x2) { + size = n; + xmin = x1; + xmax = x2; + xmaxsq = xmax*xmax; + xs = new double[n]; + ys = new double[n]; + ys1 = new double[n]; + ys2 = new double[n]; + ys3 = new double[n]; + ys4 = new double[n]; + ys5 = new double[n]; + ys6 = new double[n]; + } + virtual ~tabularFunction() { + if (xs) delete [] xs; + if (ys) delete [] ys; + if (ys1) delete [] ys1; + if (ys2) delete [] ys2; + if (ys3) delete [] ys3; + if (ys4) delete [] ys4; + if (ys5) delete [] ys5; + if (ys6) delete [] ys6; + } + void set_xrange(double x1, double x2) { + xmin = x1; + xmax = x2; + xmaxsq = xmax*xmax; + } + void set_values(int n, double x1, double x2, double * values, double epsilon) + { + int i0; + i0 = n-1; +// shrink (remove near zero points) reduces cutoff radius, and therefore computational cost +// do not shrink when x2 < 1.1 (angular function) or x2 > 20.0 (non-radial function) + if (x2 >= 1.1 && x2 <= 20.0) { + for (int i = n-1; i >= 0; i--) { + if (fabs(values[i]) > epsilon) { + i0 = i; + break; + } + } + } +// do not shrink when when list is abnormally small + if (i0 < 10/n) { + i0 = n-1; + } else if (i0 < n-1) { + values[i0] = 0.0; + i0 = i0 + 1; + values[i0] = 0.0; + } + xmin = x1; + xmax = x1 + (x2-x1)/(n -1)*i0; + xmaxsq = xmax*xmax; + n = i0+1; + resize(n); + memcpy(ys,values,n*sizeof(double)); + initialize(); + } + void value(double x, double &y, int ny, double &y1, int ny1) + { + double ps = (x - xmin) * rdx; + int ks = ps + 0.5; + if (ks > size-1) ks = size-1; + if (ks < 0 ) ks = 0; + ps = ps - ks; + if (ny) y = ((ys3[ks]*ps + ys2[ks])*ps + ys1[ks])*ps + ys[ks]; + if (ny1) y1 = (ys6[ks]*ps + ys5[ks])*ps + ys4[ks]; + } + void print_value() + { + printf("%d %f %f %f \n",size,xmin,xmax,rdx); + printf(" \n"); + for (int i = 0; i < size; i++) { + printf("%f %f \n",xs[i],ys[i]); + } + } + double get_xmin() { + return xmin; + } + double get_xmax() { + return xmax; + } + double get_xmaxsq() { + return xmaxsq; + } + double get_rdx() { + return rdx; + } + double get_vmax() { + return vmax; + } + + protected: + + void resize(int n) { + if (n != size) { + size = n; + if (xs) delete [] xs; + xs = new double[n]; + if (ys) delete [] ys; + ys = new double[n]; + if (ys1) delete [] ys1; + ys1 = new double[n]; + if (ys2) delete [] ys2; + ys2 = new double[n]; + if (ys3) delete [] ys3; + ys3 = new double[n]; + if (ys4) delete [] ys4; + ys4 = new double[n]; + if (ys5) delete [] ys5; + ys5 = new double[n]; + if (ys6) delete [] ys6; + ys6 = new double[n]; + } + } + void initialize() { + int n = size; + rdx = (xmax-xmin)/(n-1.0); + vmax = 0.0; + for (int i = 0; i < n; i++) { + if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]); + } + for (int i = 0; i < n; i++) { + xs[i] = xmin+i*rdx; + } + rdx = 1.0 / rdx; + ys1[0] = ys[1] - ys[0]; + ys1[1] = 0.5 * (ys[2] - ys[0]); + ys1[n-2] = 0.5 * (ys[n-1] - ys[n-3]); + ys1[n-1] = ys[n-1] - ys[n-2]; + for (int i = 2; i < n-2; i++) { + ys1[i]=((ys[i-2]-ys[i+2])+ 8.0*(ys[i+1]-ys[i-1]))/12.0; + } + for (int i = 0; i < n-1; i++) { + ys2[i]=3.0*(ys[i+1]-ys[i])-2.0*ys1[i]-ys1[i+1]; + ys3[i]=ys1[i]+ys1[i+1]-2.0*(ys[i+1]-ys[i]); + } + ys2[n-1]=0.0; + ys3[n-1]=0.0; + for (int i = 0; i < n; i++) { + ys4[i]=ys1[i]*rdx; + ys5[i]=2.0*ys2[i]*rdx; + ys6[i]=3.0*ys3[i]*rdx; + } + } + int size; + double xmin,xmax,xmaxsq,rdx,vmax; + double * ys, * ys1, * ys2, * ys3, * ys4, * ys5, * ys6; + double * xs; + }; + + struct PairParameters { + double cut; + double cutsq; + bool xi; // "indicator" + class tabularFunction * U; + class tabularFunction * V; + class tabularFunction * W; + class tabularFunction * P; + class tabularFunction * F; + PairParameters() { + cut = 0.0; + cutsq = 0.0; + xi = true; + U = NULL; + V = NULL; + W = NULL; + P = NULL; + F = NULL; + }; + }; + struct TripletParameters { + class tabularFunction * G; + TripletParameters() { + G = NULL; + }; + }; + + double epsilon; + bool eta; // global indicator + int nx,nr,ng; // table sizes + double maxX; + + // parameter sets + PairParameters * pairParameters; // for I-J interaction + TripletParameters * tripletParameters; // for I-J-K interaction + + int neighsize,numneighV,numneighW,numneighW1; + int *firstneighV,*firstneighW,*firstneighW1; + double *delxV,*delyV,*delzV,*drV; + double *delxW,*delyW,*delzW,*drW; + + char **elements; // names of unique elements + int **elem2param; // map: element pairs to parameters + int ***elem3param; // map: element triplets to parameters + int *map; // mapping from atom types to elements + double cutmax; // max cutoff for all elements + double cutmaxsq; + int nelements; // # of unique elements + int npair,ntriple; + int *match; + + void allocate(); + void grab(FILE *, int, double *); + + virtual void read_file(char *); + void setup(); + void write_tables(int); + + void attractive(PairParameters *, TripletParameters *, double, double, + double, double *, double *, double *, double *, double *); + + void ters_zetaterm_d(double, double *, double, double *, double, double *, + double *, double *, PairParameters *, TripletParameters *); + void costheta_d(double *, double, double *, double, + double *, double *, double *); + + // inlined functions for efficiency + + inline double vec3_dot(const double x[3], const double y[3]) const { + return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; + } + + inline void vec3_add(const double x[3], const double y[3], + double * const z) const { + z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2]; + } + + inline void vec3_scale(const double k, const double x[3], + double y[3]) const { + y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2]; + } + + inline void vec3_scaleadd(const double k, const double x[3], + const double y[3], double * const z) const { + z[0] = k*x[0]+y[0]; + z[1] = k*x[1]+y[1]; + z[2] = k*x[2]+y[2]; + } +}; + +} + +#endif +#endif diff --git a/src/USER-MISC/fix_srp.cpp b/src/USER-MISC/fix_srp.cpp index fb4bd4e5a..e5d592d30 100644 --- a/src/USER-MISC/fix_srp.cpp +++ b/src/USER-MISC/fix_srp.cpp @@ -1,609 +1,624 @@ /* ---------------------------------------------------------------------- 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 authors: Timothy Sirk (ARL), Pieter in't Veld (BASF) ------------------------------------------------------------------------- */ #include #include #include "fix_srp.h" #include "atom.h" #include "force.h" #include "domain.h" #include "comm.h" #include "memory.h" #include "error.h" #include "neighbor.h" #include "atom_vec.h" #include "modify.h" using namespace LAMMPS_NS; using namespace FixConst; /* ---------------------------------------------------------------------- */ FixSRP::FixSRP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { // settings nevery=1; peratom_freq = 1; time_integrate = 0; create_attribute = 0; comm_border = 2; // restart settings restart_global = 1; restart_peratom = 1; restart_pbc = 1; // per-atom array width 2 peratom_flag = 1; size_peratom_cols = 2; // initial allocation of atom-based array // register with Atom class array = NULL; grow_arrays(atom->nmax); // extends pack_exchange() atom->add_callback(0); atom->add_callback(1); // restart atom->add_callback(2); + // initialize to illegal values so we capture + btype = -1; + bptype = -1; + // zero for (int i = 0; i < atom->nmax; i++) for (int m = 0; m < 3; m++) array[i][m] = 0.0; - - // assume setup of fix is needed to insert particles - // is true if reading from datafile - // since a datafile cannot contain bond particles - // might be true if reading from restart - // a restart file written during the run has bond particles as per atom data - // a restart file written after the run does not have bond particles - setup = 1; } /* ---------------------------------------------------------------------- */ FixSRP::~FixSRP() { // unregister callbacks to this fix from Atom class atom->delete_callback(id,0); atom->delete_callback(id,1); atom->delete_callback(id,2); memory->destroy(array); } /* ---------------------------------------------------------------------- */ int FixSRP::setmask() { int mask = 0; mask |= PRE_FORCE; mask |= PRE_EXCHANGE; mask |= POST_RUN; return mask; } /* ---------------------------------------------------------------------- */ void FixSRP::init() { if (force->pair_match("hybrid",1) == NULL) error->all(FLERR,"Cannot use pair srp without pair_style hybrid"); - // the highest numbered atom type should be reserved for bond particles (bptype) - // set bptype, unless it will be read from restart - if(!restart_reset) bptype = atom->ntypes; - - // check if bptype is already in use - for(int i=0; i < atom->nlocal; i++) - if(atom->type[i] == bptype){ - // error if bptype is already in use - // unless starting from a rst file - // since rst can contain the bond particles as per atom data - if(!restart_reset) - error->all(FLERR,"Fix SRP requires an extra atom type"); - else - setup = 0; - } + if ((bptype < 1) || (bptype > atom->ntypes)) + error->all(FLERR,"Illegal bond particle type"); // setup neigh exclusions for diff atom types // bond particles do not interact with other types // type bptype only interacts with itself char* arg1[4]; arg1[0] = (char *) "exclude"; arg1[1] = (char *) "type"; char c0[20]; char c1[20]; - for(int z = 1; z < bptype; z++) - { - sprintf(c0, "%d", z); - arg1[2] = c0; - sprintf(c1, "%d", bptype); - arg1[3] = c1; - neighbor->modify_params(4, arg1); + for(int z = 1; z < atom->ntypes; z++) { + if(z == bptype) + continue; + sprintf(c0, "%d", z); + arg1[2] = c0; + sprintf(c1, "%d", bptype); + arg1[3] = c1; + neighbor->modify_params(4, arg1); } } /* ---------------------------------------------------------------------- insert bond particles ------------------------------------------------------------------------- */ void FixSRP::setup_pre_force(int zz) { - if(!setup) return; - - double rsqold = 0.0; - double delx, dely, delz, rmax, rsq, rsqmax; double **x = atom->x; - bigint nall = atom->nlocal + atom->nghost; double **xold; - memory->create(xold,nall,3,"fix_srp:xold"); + tagint *tag = atom->tag; + tagint *tagold; + int *type = atom->type; + int* dlist; + AtomVec *avec = atom->avec; + int **bondlist = neighbor->bondlist; + + int nlocal, nlocal_old; + nlocal = nlocal_old = atom->nlocal; + bigint nall = atom->nlocal + atom->nghost; + int nbondlist = neighbor->nbondlist; + int i,j,n; // make a copy of all coordinates and tags - // need this because create_atom overwrites ghost atoms - for(int i = 0; i < nall; i++){ + // that is consistent with the bond list as + // atom->x will be affected by creating/deleting atoms. + // also compile list of local atoms to be deleted. + + memory->create(xold,nall,3,"fix_srp:xold"); + memory->create(tagold,nall,"fix_srp:tagold"); + memory->create(dlist,nall,"fix_srp:dlist"); + + for (i = 0; i < nall; i++){ xold[i][0] = x[i][0]; xold[i][1] = x[i][1]; xold[i][2] = x[i][2]; + tagold[i]=tag[i]; + dlist[i] = (type[i] == bptype) ? 1 : 0; + for (n = 0; n < 3; n++) + array[i][n] = 0.0; } - tagint *tag = atom->tag; - tagint *tagold; - memory->create(tagold,nall,"fix_srp:tagold"); - - for(int i = 0; i < nall; i++){ - tagold[i]=tag[i]; + // delete local atoms flagged in dlist + i = 0; + int ndel = 0; + while (i < nlocal) { + if (dlist[i]) { + avec->copy(nlocal-1,i,1); + dlist[i] = dlist[nlocal-1]; + nlocal--; + ndel++; + } else i++; } - int nlocal = atom->nlocal; - int **bondlist = neighbor->bondlist; - int nbondlist = neighbor->nbondlist; + atom->nlocal = nlocal; + memory->destroy(dlist); + int nadd = 0; - int i,j; + double rsqold = 0.0; + double delx, dely, delz, rmax, rsq, rsqmax; + double xone[3]; - for (int n = 0; n < nbondlist; n++) { + for (n = 0; n < nbondlist; n++) { - // consider only the user defined bond type - // btype of zero considers all bonds - if(btype > 0 && bondlist[n][2] != btype) - continue; + // consider only the user defined bond type + // btype of zero considers all bonds + if(btype > 0 && bondlist[n][2] != btype) + continue; - i = bondlist[n][0]; - j = bondlist[n][1]; + i = bondlist[n][0]; + j = bondlist[n][1]; - // position of bond i - xone[0] = (xold[i][0] + xold[j][0])*0.5; - xone[1] = (xold[i][1] + xold[j][1])*0.5; - xone[2] = (xold[i][2] + xold[j][2])*0.5; + // position of bond i + xone[0] = (xold[i][0] + xold[j][0])*0.5; + xone[1] = (xold[i][1] + xold[j][1])*0.5; + xone[2] = (xold[i][2] + xold[j][2])*0.5; - // record longest bond - // this used to set ghost cutoff + // record longest bond + // this used to set ghost cutoff delx = xold[j][0] - xold[i][0]; dely = xold[j][1] - xold[i][1]; delz = xold[j][2] - xold[i][2]; rsq = delx*delx + dely*dely + delz*delz; if(rsq > rsqold) rsqold = rsq; - // make one particle for each bond - // i is local - // if newton bond, always make particle - // if j is local, always make particle - // if j is ghost, decide from tag - - if( force->newton_bond || j < nlocal || tagold[i] > tagold[j] ){ - atom->natoms += 1; - atom->avec->create_atom(bptype,xone); - // pack tag i/j into buffer for comm - array[atom->nlocal-1][0] = (double)tagold[i]; - array[atom->nlocal-1][1] = (double)tagold[j]; - nadd++; + // make one particle for each bond + // i is local + // if newton bond, always make particle + // if j is local, always make particle + // if j is ghost, decide from tag + + if ((force->newton_bond) || (j < nlocal_old) || (tagold[i] > tagold[j])) { + atom->natoms++; + avec->create_atom(bptype,xone); + // pack tag i/j into buffer for comm + array[atom->nlocal-1][0] = static_cast(tagold[i]); + array[atom->nlocal-1][1] = static_cast(tagold[j]); + nadd++; } } - // free temporary storage - memory->destroy(xold); - memory->destroy(tagold); - - // new total # of atoms bigint nblocal = atom->nlocal; MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); - // assign tags for new atoms, update map - if (atom->tag_enable) atom->tag_extend(); - - if (atom->map_style) { - atom->nghost = 0; - atom->map_init(); - atom->map_set(); - } + // free temporary storage + memory->destroy(xold); + memory->destroy(tagold); char str[128]; - int Nadd = 0; - MPI_Allreduce(&nadd,&Nadd,1,MPI_INT,MPI_SUM,world); + int nadd_all = 0, ndel_all = 0; + MPI_Allreduce(&ndel,&ndel_all,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&nadd,&nadd_all,1,MPI_INT,MPI_SUM,world); if(comm->me == 0){ - sprintf(str, "Inserted %d bond particles.", Nadd); + sprintf(str, "Removed/inserted %d/%d bond particles.", ndel_all,nadd_all); error->message(FLERR,str); } // check ghost comm distances // warn and change if shorter from estimate // ghost atoms must be present for bonds on edge of neighbor cutoff // extend cutghost slightly more than half of the longest bond MPI_Allreduce(&rsqold,&rsqmax,1,MPI_DOUBLE,MPI_MAX,world); rmax = sqrt(rsqmax); double cutneighmax_srp = neighbor->cutneighmax + 0.51*rmax; // find smallest cutghost double cutghostmin = comm->cutghost[0]; if (cutghostmin > comm->cutghost[1]) cutghostmin = comm->cutghost[1]; if (cutghostmin > comm->cutghost[2]) cutghostmin = comm->cutghost[2]; // reset cutghost if needed if (cutneighmax_srp > cutghostmin){ if(comm->me == 0){ sprintf(str, "Extending ghost comm cutoff. New %f, old %f.", cutneighmax_srp, cutghostmin); error->message(FLERR,str); } // cutghost updated by comm->setup comm->cutghostuser = cutneighmax_srp; } + // assign tags for new atoms, update map + atom->tag_extend(); + if (atom->map_style) { + atom->nghost = 0; + atom->map_init(); + atom->map_set(); + } + // put new particles in the box before exchange // move owned to new procs // get ghosts // build neigh lists again // if triclinic, lambda coords needed for pbc, exchange, borders if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); comm->exchange(); if (atom->sortfreq > 0) atom->sort(); comm->borders(); // back to box coords if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); domain->image_check(); domain->box_too_small_check(); modify->setup_pre_neighbor(); neighbor->build(); neighbor->ncalls = 0; // new atom counts + nlocal = atom->nlocal; nall = atom->nlocal + atom->nghost; + // zero all forces - for(int i = 0; i < nall; i++) - for(int n = 0; n < 3; n++) - atom->f[i][n] = 0.0; + + for(i = 0; i < nall; i++) + atom->f[i][0] = atom->f[i][1] = atom->f[i][2] = 0.0; // do not include bond particles in thermo output - // remove them from all groups - for(int i=0; i< nlocal; i++) - if(atom->type[i] == bptype) + // remove them from all groups. set their velocity to zero. + + for(i=0; i< nlocal; i++) + if(atom->type[i] == bptype) { atom->mask[i] = 0; + atom->v[i][0] = atom->v[i][1] = atom->v[i][2] = 0.0; + } } /* ---------------------------------------------------------------------- set position of bond particles ------------------------------------------------------------------------- */ void FixSRP::pre_exchange() { // update ghosts comm->forward_comm(); // reassign bond particle coordinates to midpoint of bonds // only need to do this before neigh rebuild double **x=atom->x; int i,j; int nlocal = atom->nlocal; for(int ii = 0; ii < nlocal; ii++){ if(atom->type[ii] != bptype) continue; - i = atom->map((int)array[ii][0]); + i = atom->map(static_cast(array[ii][0])); if(i < 0) error->all(FLERR,"Fix SRP failed to map atom"); i = domain->closest_image(ii,i); - j = atom->map((int)array[ii][1]); + j = atom->map(static_cast(array[ii][1])); if(j < 0) error->all(FLERR,"Fix SRP failed to map atom"); j = domain->closest_image(ii,j); // position of bond particle ii atom->x[ii][0] = (x[i][0] + x[j][0])*0.5; atom->x[ii][1] = (x[i][1] + x[j][1])*0.5; atom->x[ii][2] = (x[i][2] + x[j][2])*0.5; } } /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ double FixSRP::memory_usage() { double bytes = atom->nmax*2 * sizeof(double); return bytes; } /* ---------------------------------------------------------------------- allocate atom-based array ------------------------------------------------------------------------- */ void FixSRP::grow_arrays(int nmax) { memory->grow(array,nmax,2,"fix_srp:array"); array_atom = array; } /* ---------------------------------------------------------------------- copy values within local atom-based array called when move to new proc ------------------------------------------------------------------------- */ void FixSRP::copy_arrays(int i, int j, int delflag) { for (int m = 0; m < 2; m++) array[j][m] = array[i][m]; } /* ---------------------------------------------------------------------- initialize one atom's array values called when atom is created ------------------------------------------------------------------------- */ void FixSRP::set_arrays(int i) { array[i][0] = -1; array[i][1] = -1; } /* ---------------------------------------------------------------------- pack values in local atom-based array for exchange with another proc ------------------------------------------------------------------------- */ int FixSRP::pack_exchange(int i, double *buf) { for (int m = 0; m < 2; m++) buf[m] = array[i][m]; return 2; } /* ---------------------------------------------------------------------- unpack values in local atom-based array from exchange with another proc ------------------------------------------------------------------------- */ int FixSRP::unpack_exchange(int nlocal, double *buf) { for (int m = 0; m < 2; m++) array[nlocal][m] = buf[m]; return 2; } /* ---------------------------------------------------------------------- pack values for border communication at re-neighboring ------------------------------------------------------------------------- */ int FixSRP::pack_border(int n, int *list, double *buf) { // pack buf for border com int i,j; int m = 0; for (i = 0; i < n; i++) { j = list[i]; buf[m++] = array[j][0]; buf[m++] = array[j][1]; } return m; } /* ---------------------------------------------------------------------- unpack values for border communication at re-neighboring ------------------------------------------------------------------------- */ int FixSRP::unpack_border(int n, int first, double *buf) { // unpack buf into array int i,last; int m = 0; last = first + n; for (i = first; i < last; i++){ - array[i][0] = static_cast (buf[m++]); - array[i][1] = static_cast (buf[m++]); + array[i][0] = buf[m++]; + array[i][1] = buf[m++]; } return m; } /* ---------------------------------------------------------------------- remove particles after run ------------------------------------------------------------------------- */ void FixSRP::post_run() { // all bond particles are removed after each run // useful for write_data and write_restart commands // since those commands occur between runs bigint natoms_previous = atom->natoms; int nlocal = atom->nlocal; int* dlist; memory->create(dlist,nlocal,"fix_srp:dlist"); for (int i = 0; i < nlocal; i++){ if(atom->type[i] == bptype) dlist[i] = 1; else dlist[i] = 0; } // delete local atoms flagged in dlist // reset nlocal AtomVec *avec = atom->avec; int i = 0; while (i < nlocal) { if (dlist[i]) { avec->copy(nlocal-1,i,1); dlist[i] = dlist[nlocal-1]; nlocal--; } else i++; } atom->nlocal = nlocal; memory->destroy(dlist); // reset atom->natoms // reset atom->map if it exists // set nghost to 0 so old ghosts won't be mapped bigint nblocal = atom->nlocal; MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world); if (atom->map_style) { atom->nghost = 0; atom->map_init(); atom->map_set(); } // print before and after atom count bigint ndelete = natoms_previous - atom->natoms; if (comm->me == 0) { if (screen) fprintf(screen,"Deleted " BIGINT_FORMAT " atoms, new total = " BIGINT_FORMAT "\n", ndelete,atom->natoms); if (logfile) fprintf(logfile,"Deleted " BIGINT_FORMAT " atoms, new total = " BIGINT_FORMAT "\n", ndelete,atom->natoms); } // verlet calls box_too_small_check() in post_run // this check maps all bond partners // therefore need ghosts // need to convert to lambda coords before apply pbc if (domain->triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); comm->setup(); comm->exchange(); if (atom->sortfreq > 0) atom->sort(); comm->borders(); // change back to box coordinates if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); } /* ---------------------------------------------------------------------- pack values in local atom-based arrays for restart file ------------------------------------------------------------------------- */ int FixSRP::pack_restart(int i, double *buf) { int m = 0; buf[m++] = 3; buf[m++] = array[i][0]; buf[m++] = array[i][1]; return m; } /* ---------------------------------------------------------------------- unpack values from atom->extra array to restart the fix ------------------------------------------------------------------------- */ void FixSRP::unpack_restart(int nlocal, int nth) { double **extra = atom->extra; // skip to Nth set of extra values int m = 0; for (int i = 0; i < nth; i++){ - m += static_cast (extra[nlocal][m]); + m += extra[nlocal][m]; } m++; array[nlocal][0] = extra[nlocal][m++]; array[nlocal][1] = extra[nlocal][m++]; } /* ---------------------------------------------------------------------- maxsize of any atom's restart data ------------------------------------------------------------------------- */ int FixSRP::maxsize_restart() { return 3; } /* ---------------------------------------------------------------------- size of atom nlocal's restart data ------------------------------------------------------------------------- */ int FixSRP::size_restart(int nlocal) { return 3; } /* ---------------------------------------------------------------------- pack global state of Fix ------------------------------------------------------------------------- */ void FixSRP::write_restart(FILE *fp) { int n = 0; double list[3]; list[n++] = comm->cutghostuser; list[n++] = btype; list[n++] = bptype; if (comm->me == 0) { int size = n * sizeof(double); fwrite(&size,sizeof(int),1,fp); fwrite(list,sizeof(double),n,fp); } } /* ---------------------------------------------------------------------- use info from restart file to restart the Fix ------------------------------------------------------------------------- */ void FixSRP::restart(char *buf) { int n = 0; double *list = (double *) buf; comm->cutghostuser = static_cast (list[n++]); btype = static_cast (list[n++]); bptype = static_cast (list[n++]); } /* ---------------------------------------------------------------------- interface with pair class pair srp sets the bond type in this fix ------------------------------------------------------------------------- */ int FixSRP::modify_param(int narg, char **arg) { if (strcmp(arg[0],"btype") == 0) { btype = atoi(arg[1]); return 2; } + if (strcmp(arg[0],"bptype") == 0) { + bptype = atoi(arg[1]); + return 2; + } return 0; } diff --git a/src/USER-MISC/fix_srp.h b/src/USER-MISC/fix_srp.h index 5315d506d..aaeea2f03 100644 --- a/src/USER-MISC/fix_srp.h +++ b/src/USER-MISC/fix_srp.h @@ -1,72 +1,70 @@ /* -*- 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 FIX_CLASS FixStyle(SRP,FixSRP) #else #ifndef LMP_FIX_SRP_H #define LMP_FIX_SRP_H #include #include "fix.h" namespace LAMMPS_NS { class FixSRP : public Fix { public: FixSRP(class LAMMPS *, int, char **); ~FixSRP(); int setmask(); void init(); void pre_exchange(); void setup_pre_force(int); double memory_usage(); void grow_arrays(int); void copy_arrays(int, int, int); void set_arrays(int); int pack_exchange(int, double *); int unpack_exchange(int, double *); int pack_border(int, int *, double *); int unpack_border(int, int, double *); void post_run(); int pack_restart(int, double*); void unpack_restart(int, int); int maxsize_restart(); int size_restart(int); void write_restart(FILE *); void restart(char *); int modify_param(int, char **); double **array; private: - double xone[3]; int btype; int bptype; - int setup; }; } #endif #endif /* ERROR/WARNING messages: */ diff --git a/src/USER-MISC/pair_srp.cpp b/src/USER-MISC/pair_srp.cpp index 33629d845..5e254b74b 100644 --- a/src/USER-MISC/pair_srp.cpp +++ b/src/USER-MISC/pair_srp.cpp @@ -1,724 +1,744 @@ /* ---------------------------------------------------------------------- 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 authors: Timothy Sirk (ARL), Pieter in't Veld (BASF) This pair style srp command calculates a segmental repulsive force between bonds. This is useful for preventing the crossing of bonds if soft non-bonded potentials are used, such as DPD polymer chains. See the doc page for pair_style srp command for usage instructions. There is an example script for this package in examples/USER/srp. Please contact Timothy Sirk for questions (tim.sirk@us.army.mil). ------------------------------------------------------------------------- */ #include #include "pair_srp.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "memory.h" #include "error.h" #include "domain.h" #include "modify.h" #include "fix.h" #include "fix_srp.h" #include "thermo.h" #include "output.h" #include #include "citeme.h" using namespace LAMMPS_NS; #define SMALL 1.0e-10 #define BIG 1e10 #define ONETWOBIT 0x40000000 static const char cite_srp[] = "@Article{Sirk2012\n" " author = {T. Sirk and Y. Sliozberg and J. Brennan and M. Lisal and J. Andzelm},\n" " title = {An enhanced entangled polymer model for dissipative particle dynamics},\n" " journal = {J.~Chem.~Phys.},\n" " year = 2012,\n" " volume = 136,\n" " pages = {134903}\n" "}\n\n"; static int srp_instance = 0; /* ---------------------------------------------------------------------- set size of pair comms in constructor ---------------------------------------------------------------------- */ PairSRP::PairSRP(LAMMPS *lmp) : Pair(lmp) { writedata = 1; if (lmp->citeme) lmp->citeme->add(cite_srp); nextra = 1; segment = NULL; // generate unique fix-id for this pair style instance fix_id = strdup("XX_FIX_SRP"); fix_id[0] = '0' + srp_instance / 10; fix_id[1] = '0' + srp_instance % 10; ++srp_instance; // create fix SRP instance here, as it has to // be executed before all other fixes char **fixarg = new char*[3]; fixarg[0] = fix_id; fixarg[1] = (char *) "all"; fixarg[2] = (char *) "SRP"; modify->add_fix(3,fixarg); f_srp = (FixSRP *) modify->fix[modify->nfix-1]; delete [] fixarg; } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairSRP::allocate() { allocated = 1; // particles of bptype inserted by fix srp // bptype is the highest numbered atom type int n = bptype; memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); memory->create(cut, n + 1, n + 1, "pair:cut"); memory->create(a0, n + 1, n + 1, "pair:a0"); // setflag for atom types memory->create(setflag,n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; maxcount = 0; } /* ---------------------------------------------------------------------- free ------------------------------------------------------------------------- */ PairSRP::~PairSRP() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut); memory->destroy(a0); memory->destroy(segment); } // check nfix in case all fixes have already been deleted if (modify->nfix) modify->delete_fix(fix_id); free(fix_id); } /* ---------------------------------------------------------------------- compute bond-bond repulsions ------------------------------------------------------------------------- */ void PairSRP::compute(int eflag, int vflag) { // setup energy and virial if (eflag || vflag) ev_setup(eflag, vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; int i0, i1, j0, j1; int i,j,ii,jj,inum,jnum; double dijsq, dij; int *ilist,*jlist,*numneigh,**firstneigh; inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; double dx,dy,dz,ti,tj; double wd, lever0, lever1, evdwl, fpair; double fxlever0, fylever0, fzlever0, fxlever1, fylever1, fzlever1; double fx, fy, fz; evdwl = 0.0; // mapping global to local for atoms inside bond particles // exclude 1-2 neighs if requested if (neighbor->ago == 0){ remapBonds(nall); if(exclude) onetwoexclude(ilist, inum, jlist, numneigh, firstneigh); } // this pair style only used with hybrid // due to exclusions // each atom i is type bptype // each neigh j is type bptype // using midpoint distance option if(midpoint){ for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; // two atoms inside bond particle i0 = segment[i][0]; j0 = segment[i][1]; for (jj = 0; jj < jnum; jj++) { jlist = firstneigh[i]; j = jlist[jj]; // enforce 1-2 exclusions if( (sbmask(j) & exclude) ) continue; j &= NEIGHMASK; //retrieve atoms from bond particle i1 = segment[j][0]; j1 = segment[j][1]; // midpt dist bond 0 and 1 dx = 0.5*(x[i0][0] - x[i1][0] + x[j0][0] - x[j1][0]); dy = 0.5*(x[i0][1] - x[i1][1] + x[j0][1] - x[j1][1]); dz = 0.5*(x[i0][2] - x[i1][2] + x[j0][2] - x[j1][2]); dijsq = dx*dx + dy*dy + dz*dz; if (dijsq < cutsq[bptype][bptype]){ dij = sqrt(dijsq); if (dij < SMALL) continue; // dij can be 0.0 with soft potentials wd = 1.0 - dij / cut[bptype][bptype]; fpair = 0.5 * a0[bptype][bptype] * wd / dij; // 0.5 factor for lever rule // force for bond 0, beads 0,1 //force between bonds fx = fpair * dx; fy = fpair * dy; fz = fpair * dz; f[i0][0] += fx; //keep force sign for bond 0 f[i0][1] += fy; f[i0][2] += fz; f[j0][0] += fx; f[j0][1] += fy; f[j0][2] += fz; f[i1][0] -= fx; //flip force sign for bond 1 f[i1][1] -= fy; f[i1][2] -= fz; f[j1][0] -= fx; f[j1][1] -= fy; f[j1][2] -= fz; // ************************************************* // if (eflag){ evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd; } if (evflag){ ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz); ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,fpair,dx,dy,dz); } if (vflag_fdotr) virial_fdotr_compute(); } } } } else{ // using min distance option for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; i0 = segment[i][0]; j0 = segment[i][1]; for (jj = 0; jj < jnum; jj++) { jlist = firstneigh[i]; j = jlist[jj]; // enforce 1-2 exclusions if( (sbmask(j) & exclude) ) continue; j &= NEIGHMASK; i1 = segment[j][0]; j1 = segment[j][1]; getMinDist(x, dx, dy, dz, ti, tj, i0, j0, i1, j1); dijsq = dx*dx + dy*dy + dz*dz; if (dijsq < cutsq[bptype][bptype]){ dij = sqrt(dijsq); if (dij < SMALL) continue; // dij can be 0.0 with soft potentials wd = 1.0 - dij / cut[bptype][bptype]; fpair = a0[bptype][bptype] * wd / dij; // force for bond 0, beads 0,1 lever0 = 0.5 + ti; // assign force according to lever rule lever1 = 0.5 + tj; // assign force according to lever rule //force between bonds fx = fpair * dx; fy = fpair * dy; fz = fpair * dz; //decompose onto atoms fxlever0 = fx * lever0; fylever0 = fy * lever0; fzlever0 = fz * lever0; fxlever1 = fx * lever1; fylever1 = fy * lever1; fzlever1 = fz * lever1; f[i0][0] += fxlever0; //keep force sign for bond 0 f[i0][1] += fylever0; f[i0][2] += fzlever0; f[j0][0] += (fx - fxlever0); f[j0][1] += (fy - fylever0); f[j0][2] += (fz - fzlever0); f[i1][0] -= fxlever1; //flip force sign for bond 1 f[i1][1] -= fylever1; f[i1][2] -= fzlever1; f[j1][0] -= (fx - fxlever1); f[j1][1] -= (fy - fylever1); f[j1][2] -= (fz - fzlever1); // ************************************************* // if (eflag){ evdwl = 0.5 * a0[bptype][bptype] * cut[bptype][bptype] * wd * wd; } if (evflag){ ev_tally(i0,i1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz); ev_tally(j0,j1,nlocal,1,0.5*evdwl,0.0,0.5*fpair,dx,dy,dz); } if (vflag_fdotr) virial_fdotr_compute(); } } } } } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairSRP::settings(int narg, char **arg) { - if (narg < 3 || narg > 5) - error->all(FLERR,"Illegal pair_style command"); - - cut_global = force->numeric(FLERR,arg[0]); - // wildcard - if (strcmp(arg[1],"*") == 0) - btype = 0; - else { - btype = force->inumeric(FLERR,arg[1]); - if ((btype > atom->nbondtypes) || (btype <= 0)) - error->all(FLERR,"Illegal pair_style command"); - } + if (narg < 3 || narg > 7) + error->all(FLERR,"Illegal pair_style command"); + + if (atom->tag_enable == 0) + error->all(FLERR,"Pair_style srp requires atom IDs"); + + cut_global = force->numeric(FLERR,arg[0]); + // wildcard + if (strcmp(arg[1],"*") == 0) + btype = 0; + else { + btype = force->inumeric(FLERR,arg[1]); + if ((btype > atom->nbondtypes) || (btype <= 0)) + error->all(FLERR,"Illegal pair_style command"); + } - // settings - midpoint = 0; - min = 0; + // settings + midpoint = 0; + min = 0; if (strcmp(arg[2],"min") == 0) min = 1; else if (strcmp(arg[2],"mid") == 0) midpoint = 1; else - error->all(FLERR,"Illegal pair_style command"); + error->all(FLERR,"Illegal pair_style command"); int iarg = 3; // default exclude 1-2 // scaling for 1-2, etc not supported exclude = 1; + // use last atom type by default for bond particles + bptype = atom->ntypes; + while (iarg < narg) { if (strcmp(arg[iarg],"exclude") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command"); if (strcmp(arg[iarg+1],"yes") == 0) exclude = 1; if (strcmp(arg[iarg+1],"no") == 0){ if (min) error->all(FLERR,"Illegal exclude option in pair srp command"); exclude = 0; } iarg += 2; + } else if (strcmp(arg[iarg],"bptype") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal pair srp command"); + bptype = force->inumeric(FLERR,arg[iarg+1]); + if ((bptype < 1) || (bptype > atom->ntypes)) + error->all(FLERR,"Illegal bond particle type for srp"); + iarg += 2; } else error->all(FLERR,"Illegal pair srp command"); } - // highest atom type is saved for bond particles - bptype = atom->ntypes; - // reset cutoffs if explicitly set if (allocated) { int i,j; for (i = 1; i <= bptype; i++) for (j = i+1; j <= bptype; j++) if (setflag[i][j]) cut[i][j] = cut_global; } } /* ---------------------------------------------------------------------- set coeffs ------------------------------------------------------------------------- */ void PairSRP::coeff(int narg, char **arg) { if (narg < 3 || narg > 4) error->all(FLERR,"PairSRP: Incorrect args for pair coeff"); if (!allocated) allocate(); // set ij bond-bond cutoffs int ilo, ihi, jlo, jhi; force->bounds(arg[0], bptype, ilo, ihi); force->bounds(arg[1], bptype, jlo, jhi); double a0_one = force->numeric(FLERR,arg[2]); double cut_one = cut_global; if (narg == 4) cut_one = force->numeric(FLERR,arg[3]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { a0[i][j] = a0_one; cut[i][j] = cut_one; cutsq[i][j] = cut_one * cut_one; setflag[i][j] = 1; count++; } } if (count == 0) error->warning(FLERR,"PairSRP: No pair coefficients were set"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairSRP::init_style() { if (!force->newton_pair) error->all(FLERR,"PairSRP: Pair srp requires newton pair on"); // verify that fix SRP is still defined and has not been changed. int ifix = modify->find_fix(fix_id); if (f_srp != (FixSRP *)modify->fix[ifix]) error->all(FLERR,"Fix SRP has been changed unexpectedly"); - // set bond type in fix srp + if (comm->me == 0) { + if (screen) fprintf(screen,"Using type %d for bond particles\n",bptype); + if (logfile) fprintf(logfile,"Using type %d for bond particles\n",bptype); + } + + // set bond and bond particle types in fix srp // bonds of this type will be represented by bond particles + // if bond type is 0, then all bonds have bond particles // btype = bond type - // bptype = bond particle type char c0[20]; char* arg0[2]; sprintf(c0, "%d", btype); arg0[0] = (char *) "btype"; arg0[1] = c0; f_srp->modify_params(2, arg0); + // bptype = bond particle type + sprintf(c0, "%d", bptype); + arg0[0] = (char *) "bptype"; + arg0[1] = c0; + f_srp->modify_params(2, arg0); + // bond particles do not contribute to energy or virial // bond particles do not belong to group all // but thermo normalization is by nall // therefore should turn off normalization int me; MPI_Comm_rank(world,&me); char *arg1[2]; arg1[0] = (char *) "norm"; arg1[1] = (char *) "no"; output->thermo->modify_params(2, arg1); if (me == 0) error->message(FLERR,"Thermo normalization turned off by pair srp"); neighbor->request(this,instance_me); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairSRP::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"PairSRP: All pair coeffs are not set"); cut[j][i] = cut[i][j]; a0[j][i] = a0[i][j]; return cut[i][j]; } /* ---------------------------------------------------------------------- find min distance for bonds i0/j0 and i1/j1 ------------------------------------------------------------------------- */ inline void PairSRP::getMinDist(double** &x, double &dx, double &dy, double &dz, double &ti, double &tj, int &i0, int &j0, int &i1, int &j1) { // move these outside the loop double diffx0, diffy0, diffz0, diffx1, diffy1, diffz1, dPx, dPy, dPz, RiRi, RiRj, RjRj; double denom, termx0, termy0, termz0, num0, termx1, termy1, termz1, num1; // compute midpt dist from 1st atom, 1st bond diffx0 = x[j0][0] - x[i0][0]; // x,y,z from bond 0 diffy0 = x[j0][1] - x[i0][1]; diffz0 = x[j0][2] - x[i0][2]; // compute midpt dist from 1st atom, 2nd bond diffx1 = x[j1][0] - x[i1][0]; diffy1 = x[j1][1] - x[i1][1]; diffz1 = x[j1][2] - x[i1][2]; // midpoint distance dPx = 0.5*(diffx0-diffx1) + x[i0][0]-x[i1][0]; dPy = 0.5*(diffy0-diffy1) + x[i0][1]-x[i1][1]; dPz = 0.5*(diffz0-diffz1) + x[i0][2]-x[i1][2]; // Ri^2 Rj^2 RiRi = diffx0*diffx0 + diffy0*diffy0 + diffz0*diffz0; RiRj = diffx0*diffx1 + diffy0*diffy1 + diffz0*diffz1; RjRj = diffx1*diffx1 + diffy1*diffy1 + diffz1*diffz1; denom = RiRj*RiRj - RiRi*RjRj; // handle case of parallel lines // reduce to midpt distance if (fabs(denom) < SMALL){ if(denom < 0) denom = -BIG; else denom = BIG; } // calc ti termx0 = RiRj*diffx1 - RjRj*diffx0; termy0 = RiRj*diffy1 - RjRj*diffy0; termz0 = RiRj*diffz1 - RjRj*diffz0; num0 = dPx*termx0 + dPy*termy0 + dPz*termz0; ti = num0 / denom; if (ti > 0.5) ti = 0.5; if (ti < -0.5) ti = -0.5; // calc tj termx1 = RiRj*diffx0 - RiRi*diffx1; termy1 = RiRj*diffy0 - RiRi*diffy1; termz1 = RiRj*diffz0 - RiRi*diffz1; num1 = dPx*termx1 + dPy*termy1 + dPz*termz1; tj = -num1/ denom; if (tj > 0.5) tj = 0.5; if (tj < -0.5) tj = -0.5; // min dist dx = dPx - ti*diffx0 + tj*diffx1; dy = dPy - ti*diffy0 + tj*diffy1; dz = dPz - ti*diffz0 + tj*diffz1; } /* -------------------------------------------------------- map global id of atoms in stored by each bond particle ------------------------------------------------------- */ inline void PairSRP::remapBonds(int &nall) { if(nall > maxcount){ memory->grow(segment, nall, 2, "pair:segment"); maxcount = nall; } // loop over all bond particles // each bond paricle holds two bond atoms // map global ids of bond atoms to local ids // might not be able to map both bond atoms of j, if j is outside neighcut // these are not on neighlist, so are not used int tmp; srp = f_srp->array_atom; for (int i = 0; i < nall; i++) { if(atom->type[i] == bptype){ // tmp is local id // tmp == -1 is ok tmp = atom->map((int)srp[i][0]); segment[i][0] = domain->closest_image(i,tmp); // repeat with other id tmp = atom->map((int)srp[i][1]); segment[i][1] = domain->closest_image(i,tmp); } } } /* -------------------------------------------------------- add exclusions for 1-2 neighs, if requested more complex exclusions or scaling probably not needed ------------------------------------------------------- */ inline void PairSRP::onetwoexclude(int* &ilist, int &inum, int* &jlist, int* &numneigh, int** &firstneigh) { int i0, i1, j0, j1; int i,j,ii,jj,jnum; // encode neighs with exclusions // only need 1-2 info for normal uses of srp // add 1-3, etc later if ever needed for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; // two atoms inside bond particle i0 = segment[i][0]; j0 = segment[i][1]; for (jj = 0; jj < jnum; jj++) { jlist = firstneigh[i]; j = jlist[jj]; j &= NEIGHMASK; //two atoms inside bond particle i1 = segment[j][0]; j1 = segment[j][1]; // check for a 1-2 neigh if(i0 == i1 || i0 == j1 || i1 == j0 || j0 == j1){ j |= ONETWOBIT; jlist[jj] = j; } } } } /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairSRP::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g\n",i,a0[i][i]); } /* ---------------------------------------------------------------------- proc 0 writes all pairs to data file ------------------------------------------------------------------------- */ void PairSRP::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) fprintf(fp,"%d %d %g %g\n",i,j,a0[i][j],cut[i][j]); } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairSRP::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&a0[i][j],sizeof(double),1,fp); fwrite(&cut[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairSRP::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { printf(" i %d j %d \n",i,j); fread(&a0[i][j],sizeof(double),1,fp); fread(&cut[i][j],sizeof(double),1,fp); } MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairSRP::write_restart_settings(FILE *fp) { fwrite(&cut_global,sizeof(double),1,fp); fwrite(&bptype,sizeof(int),1,fp); fwrite(&btype,sizeof(int),1,fp); fwrite(&min,sizeof(int),1,fp); fwrite(&midpoint,sizeof(int),1,fp); fwrite(&exclude,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairSRP::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&cut_global,sizeof(double),1,fp); fread(&bptype,sizeof(int),1,fp); fread(&btype,sizeof(int),1,fp); fread(&min,sizeof(int),1,fp); fread(&midpoint,sizeof(int),1,fp); fread(&exclude,sizeof(int),1,fp); } MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); }