diff --git a/src/MANYBODY/fix_qeq_comb.cpp b/src/MANYBODY/fix_qeq_comb.cpp
index e01a8cefb..b734c5a1c 100644
--- a/src/MANYBODY/fix_qeq_comb.cpp
+++ b/src/MANYBODY/fix_qeq_comb.cpp
@@ -1,255 +1,255 @@
-/* ----------------------------------------------------------------------
- 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: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
-------------------------------------------------------------------------- */
-
-#include "mpi.h"
-#include "math.h"
-#include "stdlib.h"
-#include "string.h"
-#include "fix_qeq_comb.h"
-#include "lmptype.h"
-#include "atom.h"
-#include "force.h"
-#include "group.h"
-#include "respa.h"
-#include "pair_comb.h"
-#include "update.h"
-#include "memory.h"
-#include "error.h"
-
-using namespace LAMMPS_NS;
-
-#define MIN(A,B) ((A) < (B)) ? (A) : (B)
-#define MAX(A,B) ((A) > (B)) ? (A) : (B)
-
-/* ---------------------------------------------------------------------- */
-
-FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
-{
- if (narg < 5) error->all("Illegal fix qeq/comb command");
-
- peratom_flag = 1;
- size_peratom_cols = 0;
- peratom_freq = 1;
-
- nevery = force->inumeric(arg[3]);
- precision = force->numeric(arg[4]);
-
- if (nevery <= 0 || precision <= 0.0)
- error->all("Illegal fix qeq/comb command");
-
- MPI_Comm_rank(world,&me);
-
- // optional args
-
- fp = NULL;
-
- int iarg = 5;
- while (iarg < narg) {
- if (strcmp(arg[iarg],"file") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix qeq/comb command");
- if (me == 0) {
- fp = fopen(arg[iarg+1],"w");
- if (fp == NULL) {
- char str[128];
- sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]);
- error->one(str);
- }
- }
- iarg += 2;
- } else error->all("Illegal fix qeq/comb command");
- }
-
- nmax = atom->nmax;
- qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
- q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
- q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
- vector_atom = qf;
-
- // zero the vector since dump may access it on timestep 0
- // zero the vector since a variable may access it before first run
-
- int nlocal = atom->nlocal;
- for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-FixQEQComb::~FixQEQComb()
-{
- if (me == 0 && fp) fclose(fp);
- memory->sfree(qf);
- memory->sfree(q1);
- memory->sfree(q2);
-}
-
-/* ---------------------------------------------------------------------- */
-
-int FixQEQComb::setmask()
-{
- int mask = 0;
- mask |= POST_FORCE;
- mask |= POST_FORCE_RESPA;
- return mask;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixQEQComb::init()
-{
- if (!atom->q_flag)
- error->all("Fix qeq/comb requires atom attribute q");
-
- comb = (PairComb *) force->pair_match("comb",1);
- if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb");
-
- if (strcmp(update->integrate_style,"respa") == 0)
- nlevels_respa = ((Respa *) update->integrate)->nlevels;
-
- ngroup = group->count(igroup);
- if (ngroup == 0) error->all("Fix qeq/comb group has no atoms");
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixQEQComb::setup(int vflag)
-{
- firstflag = 1;
- if (strcmp(update->integrate_style,"verlet") == 0)
- post_force(vflag);
- else {
- ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
- post_force_respa(vflag,nlevels_respa-1,0);
- ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
- }
- firstflag = 0;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixQEQComb::post_force(int vflag)
-{
- int i,iloop,loopmax;
- double heatpq,qmass,dtq,dtq2;
- double enegchkall,enegmaxall;
-
- if (update->ntimestep % nevery) return;
-
- // reallocate work arrays if necessary
- // qf = charge force
- // q1 = charge displacement
- // q2 = tmp storage of charge force for next iteration
-
- if (atom->nmax > nmax) {
- memory->sfree(qf);
- memory->sfree(q1);
- memory->sfree(q2);
- nmax = atom->nmax;
- qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
- q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
- q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
- vector_atom = qf;
- }
-
- // more loops for first-time charge equilibrium
-
- iloop = 0;
- if (firstflag) loopmax = 1000;
- else loopmax = 500;
-
- // charge-equilibration loop
-
- if (me == 0 && fp)
- fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
- update->ntimestep);
-
- heatpq = 0.01;
- qmass = 0.06;
- dtq = 0.040;
- dtq2 = 0.5*dtq*dtq/qmass;
-
- double enegchk = 0.0;
- double enegtot = 0.0;
- double enegmax = 0.0;
-
- double *q = atom->q;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
-
- for (i = 0; i < nlocal; i++)
- q1[i] = q2[i] = qf[i] = 0.0;
-
- for (iloop = 0; iloop < loopmax; iloop ++ ) {
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit) {
- q1[i] = qf[i]*dtq2 - heatpq*q1[i];
- q[i] += q1[i];
- }
-
- enegtot = comb->yasu_char(qf,igroup);
- enegtot /= ngroup;
- enegchk = enegmax = 0.0;
-
- for (i = 0; i < nlocal ; i++)
- if (mask[i] & groupbit) {
- q2[i] = enegtot-qf[i];
- enegmax = MAX(enegmax,fabs(q2[i]));
- enegchk += fabs(q2[i]);
- qf[i] = q2[i];
- }
-
- MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
- enegchk = enegchkall/ngroup;
- MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
- enegmax = enegmaxall;
-
- if (enegchk <= precision && enegmax <= 10.0*precision) break;
-
- if (me == 0 && fp)
- fprintf(fp," iteration: %d, enegtot %.6g, "
- "enegmax %.6g, fq deviation: %.6g\n",
- iloop,enegtot,enegmax,enegchk);
-
- for (i = 0; i < nlocal; i++)
- if (mask[i] & groupbit)
- q1[i] += qf[i]*dtq2 - heatpq*q1[i];
- }
-
- if (me == 0 && fp) {
- if (iloop == loopmax)
- fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
- else
- fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
- iloop,enegchk);
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
-{
- if (ilevel == nlevels_respa-1) post_force(vflag);
-}
-
-/* ----------------------------------------------------------------------
- memory usage of local atom-based arrays
-------------------------------------------------------------------------- */
-
-double FixQEQComb::memory_usage()
-{
- double bytes = atom->nmax*3 * sizeof(double);
- return bytes;
-}
+/* ----------------------------------------------------------------------
+ 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: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "math.h"
+#include "stdlib.h"
+#include "string.h"
+#include "fix_qeq_comb.h"
+#include "lmptype.h"
+#include "atom.h"
+#include "force.h"
+#include "group.h"
+#include "respa.h"
+#include "pair_comb.h"
+#include "update.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MIN(A,B) ((A) < (B)) ? (A) : (B)
+#define MAX(A,B) ((A) > (B)) ? (A) : (B)
+
+/* ---------------------------------------------------------------------- */
+
+FixQEQComb::FixQEQComb(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
+{
+ if (narg < 5) error->all("Illegal fix qeq/comb command");
+
+ peratom_flag = 1;
+ size_peratom_cols = 0;
+ peratom_freq = 1;
+
+ nevery = force->inumeric(arg[3]);
+ precision = force->numeric(arg[4]);
+
+ if (nevery <= 0 || precision <= 0.0)
+ error->all("Illegal fix qeq/comb command");
+
+ MPI_Comm_rank(world,&me);
+
+ // optional args
+
+ fp = NULL;
+
+ int iarg = 5;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"file") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix qeq/comb command");
+ if (me == 0) {
+ fp = fopen(arg[iarg+1],"w");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open fix qeq/comb file %s",arg[iarg+1]);
+ error->one(str);
+ }
+ }
+ iarg += 2;
+ } else error->all("Illegal fix qeq/comb command");
+ }
+
+ nmax = atom->nmax;
+ qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
+ q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
+ q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
+ vector_atom = qf;
+
+ // zero the vector since dump may access it on timestep 0
+ // zero the vector since a variable may access it before first run
+
+ int nlocal = atom->nlocal;
+ for (int i = 0; i < nlocal; i++) qf[i] = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixQEQComb::~FixQEQComb()
+{
+ if (me == 0 && fp) fclose(fp);
+ memory->sfree(qf);
+ memory->sfree(q1);
+ memory->sfree(q2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixQEQComb::setmask()
+{
+ int mask = 0;
+ mask |= POST_FORCE;
+ mask |= POST_FORCE_RESPA;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::init()
+{
+ if (!atom->q_flag)
+ error->all("Fix qeq/comb requires atom attribute q");
+
+ comb = (PairComb *) force->pair_match("comb",1);
+ if (comb == NULL) error->all("Must use pair_style comb with fix qeq/comb");
+
+ if (strcmp(update->integrate_style,"respa") == 0)
+ nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
+ ngroup = group->count(igroup);
+ if (ngroup == 0) error->all("Fix qeq/comb group has no atoms");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::setup(int vflag)
+{
+ firstflag = 1;
+ if (strcmp(update->integrate_style,"verlet") == 0)
+ post_force(vflag);
+ else {
+ ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
+ post_force_respa(vflag,nlevels_respa-1,0);
+ ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
+ }
+ firstflag = 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::post_force(int vflag)
+{
+ int i,iloop,loopmax;
+ double heatpq,qmass,dtq,dtq2;
+ double enegchkall,enegmaxall;
+
+ if (update->ntimestep % nevery) return;
+
+ // reallocate work arrays if necessary
+ // qf = charge force
+ // q1 = charge displacement
+ // q2 = tmp storage of charge force for next iteration
+
+ if (atom->nmax > nmax) {
+ memory->sfree(qf);
+ memory->sfree(q1);
+ memory->sfree(q2);
+ nmax = atom->nmax;
+ qf = (double *) memory->smalloc(nmax*sizeof(double),"qeq:qf");
+ q1 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q1");
+ q2 = (double *) memory->smalloc(nmax*sizeof(double),"qeq:q2");
+ vector_atom = qf;
+ }
+
+ // more loops for first-time charge equilibrium
+
+ iloop = 0;
+ if (firstflag) loopmax = 5000;
+ else loopmax = 2000;
+
+ // charge-equilibration loop
+
+ if (me == 0 && fp)
+ fprintf(fp,"Charge equilibration on step " BIGINT_FORMAT "\n",
+ update->ntimestep);
+
+ heatpq = 0.05;
+ qmass = 0.000548580;
+ dtq = 0.0006;
+ dtq2 = 0.5*dtq*dtq/qmass;
+
+ double enegchk = 0.0;
+ double enegtot = 0.0;
+ double enegmax = 0.0;
+
+ double *q = atom->q;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+
+ for (i = 0; i < nlocal; i++)
+ q1[i] = q2[i] = qf[i] = 0.0;
+
+ for (iloop = 0; iloop < loopmax; iloop ++ ) {
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ q1[i] += qf[i]*dtq2 - heatpq*q1[i];
+ q[i] += q1[i];
+ }
+
+ enegtot = comb->yasu_char(qf,igroup);
+ enegtot /= ngroup;
+ enegchk = enegmax = 0.0;
+
+ for (i = 0; i < nlocal ; i++)
+ if (mask[i] & groupbit) {
+ q2[i] = enegtot-qf[i];
+ enegmax = MAX(enegmax,fabs(q2[i]));
+ enegchk += fabs(q2[i]);
+ qf[i] = q2[i];
+ }
+
+ MPI_Allreduce(&enegchk,&enegchkall,1,MPI_DOUBLE,MPI_SUM,world);
+ enegchk = enegchkall/ngroup;
+ MPI_Allreduce(&enegmax,&enegmaxall,1,MPI_DOUBLE,MPI_MAX,world);
+ enegmax = enegmaxall;
+
+ if (enegchk <= precision && enegmax <= 100.0*precision) break;
+
+ if (me == 0 && fp)
+ fprintf(fp," iteration: %d, enegtot %.6g, "
+ "enegmax %.6g, fq deviation: %.6g\n",
+ iloop,enegtot,enegmax,enegchk);
+
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit)
+ q1[i] += qf[i]*dtq2 - heatpq*q1[i];
+ }
+
+ if (me == 0 && fp) {
+ if (iloop == loopmax)
+ fprintf(fp,"Charges did not converge in %d iterations\n",iloop);
+ else
+ fprintf(fp,"Charges converged in %d iterations to %.10f tolerance\n",
+ iloop,enegchk);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixQEQComb::post_force_respa(int vflag, int ilevel, int iloop)
+{
+ if (ilevel == nlevels_respa-1) post_force(vflag);
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double FixQEQComb::memory_usage()
+{
+ double bytes = atom->nmax*3 * sizeof(double);
+ return bytes;
+}
diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp
index 96f8159aa..c0a7bbfeb 100644
--- a/src/MANYBODY/pair_comb.cpp
+++ b/src/MANYBODY/pair_comb.cpp
@@ -1,2032 +1,2031 @@
-/* ----------------------------------------------------------------------
- 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: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
- LAMMPS implementation of the Charge-optimized many-body (COMB) potential
- based on the HELL MD program (Prof Simon Phillpot, UF, sphil@mse.ufl.edu)
- and Aidan Thompson's Tersoff code in LAMMPS
-------------------------------------------------------------------------- */
-
-#include "math.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "pair_comb.h"
-#include "atom.h"
-#include "comm.h"
-#include "force.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
-#include "memory.h"
-#include "error.h"
-#include "group.h"
-
-using namespace LAMMPS_NS;
-
-#define MAXLINE 1024
-#define DELTA 4
-
-/* ---------------------------------------------------------------------- */
-
-PairComb::PairComb(LAMMPS *lmp) : Pair(lmp)
-{
- single_enable = 0;
- one_coeff = 1;
-
- PI = 4.0*atan(1.0);
- PI2 = 2.0*atan(1.0);
- PI4 = atan(1.0);
- PIsq = sqrt(PI);
-
- nmax = 0;
- NCo = NULL;
-
- nelements = 0;
- elements = NULL;
- nparams = 0;
- maxparam = 0;
- params = NULL;
- elem2param = NULL;
-
- intype = NULL;
- fafb = NULL;
- dfafb = NULL;
- ddfafb = NULL;
- phin = NULL;
- dphin = NULL;
- erpaw = NULL;
-
- // set comm size needed by this Pair
-
- comm_forward = 1;
- comm_reverse = 1;
-}
-
-/* ----------------------------------------------------------------------
- check if allocated, since class can be destructed when incomplete
-------------------------------------------------------------------------- */
-
-PairComb::~PairComb()
-{
- memory->sfree(NCo);
-
- if (elements)
- for (int i = 0; i < nelements; i++) delete [] elements[i];
- delete [] elements;
- memory->sfree(params);
- memory->destroy_3d_int_array(elem2param);
-
- memory->destroy_2d_int_array(intype);
- memory->destroy_2d_double_array(fafb);
- memory->destroy_2d_double_array(dfafb);
- memory->destroy_2d_double_array(ddfafb);
- memory->destroy_2d_double_array(phin);
- memory->destroy_2d_double_array(dphin);
- memory->destroy_2d_double_array(erpaw);
-
- if (allocated) {
- memory->destroy_2d_int_array(setflag);
- memory->destroy_2d_double_array(cutsq);
- delete [] map;
- delete [] esm;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::compute(int eflag, int vflag)
-{
- int i,j,k,ii,jj,kk,inum,jnum,iparam_i;
- int itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk;
- double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
- double rsq,rsq1,rsq2;
- double delr1[3],delr2[3],fi[3],fj[3],fk[3];
- double zeta_ij,prefactor;
- double fqi,fqij;
- int *ilist,*jlist,*numneigh,**firstneigh;
- int mr1,mr2,mr3;
- int rsc,inty;
- double elp_ij,filp[3],fjlp[3],fklp[3];
- double iq,jq;
- double yaself;
- double potal,fac11,fac11e;
- double vionij,fvionij,sr1,sr2,sr3,Eov,Fov;
-
- evdwl = ecoul = 0.0;
- if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = vflag_fdotr = vflag_atom = 0;
-
- // grow coordination array if necessary
-
- if (atom->nmax > nmax) {
- memory->sfree(NCo);
- nmax = atom->nmax;
- NCo = (int *) memory->smalloc(nmax*sizeof(double),"pair:NCo");
- }
-
- double **x = atom->x;
- double **f = atom->f;
- double *q = atom->q;
- int *tag = atom->tag;
- int *type = atom->type;
- int nlocal = atom->nlocal;
- int newton_pair = force->newton_pair;
- int ntype = atom->ntypes;
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- yaself = vionij = fvionij = Eov = Fov = 0.0;
-
- // self energy correction term: potal
-
- potal_calc(potal,fac11,fac11e);
-
- // 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];
- iq = q[i];
- NCo[i] = 0;
- iparam_i = elem2param[itype][itype][itype];
-
- // self energy, only on i atom
-
- yaself = self(¶ms[iparam_i],iq,potal);
-
- if (evflag) ev_tally(i,i,nlocal,0,yaself,0.0,0.0,0.0,0.0,0.0);
-
- // two-body interactions (long and short repulsive)
-
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- 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;
- }
-
- // Qj calculates 2-body Coulombic
-
- jtype = map[type[j]];
- jq = q[j];
-
- 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[itype][jtype][jtype];
-
- // long range q-dependent
-
- if (rsq > params[iparam_ij].lcutsq) continue;
-
- inty = intype[itype][jtype];
-
- // polynomial three-point interpolation
-
- tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3, itype);
-
- // 1/r energy and forces
-
- direct(inty,mr1,mr2,mr3,rsq,sr1,sr2,sr3,iq,jq,
- potal,fac11,fac11e,vionij,fvionij);
-
- // field correction to self energy
-
- field(¶ms[iparam_ij],rsq,iq,jq,vionij,fvionij);
-
- // polarization field
- // sums up long range forces
-
- f[i][0] += delx*fvionij;
- f[i][1] += dely*fvionij;
- f[i][2] += delz*fvionij;
- f[j][0] -= delx*fvionij;
- f[j][1] -= dely*fvionij;
- f[j][2] -= delz*fvionij;
-
- if (evflag)
- ev_tally(i,j,nlocal,newton_pair,0.0,vionij,fvionij,delx,dely,delz);
-
- // short range q-independent
-
- if (rsq > params[iparam_ij].cutsq) continue;
-
- repulsive(¶ms[iparam_ij],rsq,fpair,eflag,evdwl,iq,jq);
-
- // repulsion is pure two-body, sums up pair repulsive forces
-
- 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);
- }
-
- // accumulate coordination number information
-
- if (cor_flag) {
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- jtype = map[type[j]];
- iparam_ij = elem2param[itype][jtype][jtype];
-
- if(params[iparam_ij].hfocor > 0.0 ) {
- delr1[0] = x[j][0] - xtmp;
- delr1[1] = x[j][1] - ytmp;
- delr1[2] = x[j][2] - ztmp;
- rsq1 = vec3_dot(delr1,delr1);
-
- if (rsq1 > params[iparam_ij].cutsq) continue;
- NCo[i] += 1;
- }
- }
- }
-
- // three-body interactions
- // skip immediately if I-J is not within cutoff
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- jtype = map[type[j]];
- iparam_ij = elem2param[itype][jtype][jtype];
-
- // this Qj for q-dependent BSi
-
- jq = q[j];
-
- delr1[0] = x[j][0] - xtmp;
- delr1[1] = x[j][1] - ytmp;
- delr1[2] = x[j][2] - ztmp;
- rsq1 = vec3_dot(delr1,delr1);
-
- if (rsq1 > params[iparam_ij].cutsq) continue;
-
- // accumulate bondorder zeta for each i-j interaction via loop over k
-
- zeta_ij = 0.0;
- cuo_flag1 = 0; cuo_flag2 = 0;
-
- for (kk = 0; kk < jnum; kk++) {
- if (jj == kk) continue;
- k = jlist[kk];
- ktype = map[type[k]];
- iparam_ijk = elem2param[itype][jtype][ktype];
-
- delr2[0] = x[k][0] - xtmp;
- delr2[1] = x[k][1] - ytmp;
- delr2[2] = x[k][2] - ztmp;
- rsq2 = vec3_dot(delr2,delr2);
-
- if (rsq2 > params[iparam_ijk].cutsq) continue;
-
- zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
-
- if (params[iparam_ijk].hfocor == -2.0) cuo_flag1 = 1;
- if (params[iparam_ijk].hfocor == -1.0) cuo_flag2 = 1;
- }
-
- if (cuo_flag1 && cuo_flag2) cuo_flag = 1;
- else cuo_flag = 0;
-
- force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fpair,
- prefactor,eflag,evdwl,iq,jq);
-
- // over-coordination correction for HfO2
-
- if (cor_flag && NCo[i] != 0)
- Over_cor(¶ms[iparam_ij],rsq1,NCo[i],Eov, Fov);
- evdwl += Eov;
- fpair += Fov;
-
- 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 (3-body forces)
-
- for (kk = 0; kk < jnum; kk++) {
- if (jj == kk) continue;
- k = jlist[kk];
- ktype = map[type[k]];
- iparam_ijk = elem2param[itype][jtype][ktype];
-
- delr2[0] = x[k][0] - xtmp;
- delr2[1] = x[k][1] - ytmp;
- delr2[2] = x[k][2] - ztmp;
- rsq2 = vec3_dot(delr2,delr2);
- if (rsq2 > params[iparam_ijk].cutsq) continue;
-
- for (rsc = 0; rsc < 3; rsc++)
- fi[rsc] = fj[rsc] = fk[rsc] = 0.0;
-
- attractive(¶ms[iparam_ijk],prefactor,
- rsq1,rsq2,delr1,delr2,fi,fj,fk);
-
- // 3-body LP and BB correction and forces
-
- elp_ij = elp(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
- flp(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2,filp,fjlp,fklp);
-
- for (rsc = 0; rsc < 3; rsc++) {
- fi[rsc] += filp[rsc];
- fj[rsc] += fjlp[rsc];
- fk[rsc] += fklp[rsc];
- }
-
- for (rsc = 0; rsc < 3; rsc++) {
- f[i][rsc] += fi[rsc];
- f[j][rsc] += fj[rsc];
- f[k][rsc] += fk[rsc];
- }
-
- if (evflag)
- ev_tally(i,j,nlocal,newton_pair,elp_ij,0.0,0.0,0.0,0.0,0.0);
- if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);
-
- }
- }
-
- if (cuo_flag) params[iparam_i].cutsq *= 0.65;
- }
-
- cuo_flag = 0;
-
- if (vflag_fdotr) virial_compute();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::allocate()
-{
- allocated = 1;
- int n = atom->ntypes;
-
- setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
- cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
-
- map = new int[n+1];
- esm = new double[n];
-}
-
-/* ----------------------------------------------------------------------
- global settings
-------------------------------------------------------------------------- */
-
-void PairComb::settings(int narg, char **arg)
-{
- if (narg > 0) error->all("Illegal pair_style command");
-}
-
-/* ----------------------------------------------------------------------
- set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairComb::coeff(int narg, char **arg)
-{
- int i,j,n;
- double r;
-
- if (!allocated) allocate();
-
- if (narg != 3 + atom->ntypes)
- error->all("Incorrect args for pair coefficients");
-
- // insure I,J args are * *
-
- if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
- error->all("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();
-
- n = atom->ntypes;
-
- // generate streitz-mintmire direct 1/r energy look-up table
-
- if (comm->me == 0 && screen) fprintf(screen,"Pair COMB:\n");
- if (comm->me == 0 && screen)
- fprintf(screen," generating Coulomb integral lookup table ...\n");
- sm_table();
-
- if (cor_flag && comm->me == 0 && screen)
- fprintf(screen," will apply over-coordination correction ...\n");
- if (!cor_flag&& comm->me == 0 && screen)
- fprintf(screen," will not apply over-coordination correction ...\n");
-
- // clear setflag since coeff() called once with I,J = * *
-
- 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("Incorrect args for pair coefficients");
-}
-
-/* ----------------------------------------------------------------------
- init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairComb::init_style()
-{
- if (atom->tag_enable == 0)
- error->all("Pair style COMB requires atom IDs");
- if (force->newton_pair == 0)
- error->all("Pair style COMB requires newton pair on");
- if (!atom->q_flag)
- error->all("Pair style COMB requires atom attribute q");
-
- // ptr to QEQ fix
-
- //for (i = 0; i < modify->nfix; i++)
- // if (strcmp(modify->fix[i]->style,"qeq") == 0) break;
- //if (i < modify->nfix) fixqeq = (FixQEQ *) modify->fix[i];
- //else fixqeq = NULL;
-
- // 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 PairComb::init_one(int i, int j)
-{
- if (setflag[i][j] == 0) error->all("All pair coeffs are not set");
- return cutmax;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::read_file(char *file)
-{
- int params_per_line = 49;
- char **words = new char*[params_per_line+1];
-
- if (params) delete [] params;
- params = NULL;
- nparams = 0;
-
- // open file on proc 0
-
- FILE *fp;
- if (comm->me == 0) {
- fp = fopen(file,"r");
- if (fp == NULL) {
- char str[128];
- sprintf(str,"Cannot open COMB potential file %s",file);
- error->one(str);
- }
- }
-
- // read each line out of file, skipping blank lines or leading '#'
- // store line of params if all 3 element tags are in element list
-
- int n,nwords,ielement,jelement,kelement;
- char line[MAXLINE],*ptr;
- int eof = 0;
-
- while (1) {
- if (comm->me == 0) {
- ptr = fgets(line,MAXLINE,fp);
- if (ptr == NULL) {
- eof = 1;
- fclose(fp);
- } else n = strlen(line) + 1;
- }
- MPI_Bcast(&eof,1,MPI_INT,0,world);
- if (eof) break;
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
-
- // strip comment, skip line if blank
-
- if (ptr = strchr(line,'#')) *ptr = '\0';
- nwords = atom->count_words(line);
- if (nwords == 0) continue;
-
- // concatenate additional lines until have params_per_line words
-
- while (nwords < params_per_line) {
- n = strlen(line);
- if (comm->me == 0) {
- ptr = fgets(&line[n],MAXLINE-n,fp);
- if (ptr == NULL) {
- eof = 1;
- fclose(fp);
- } else n = strlen(line) + 1;
- }
- MPI_Bcast(&eof,1,MPI_INT,0,world);
- if (eof) break;
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- if (ptr = strchr(line,'#')) *ptr = '\0';
- nwords = atom->count_words(line);
- }
-
- if (nwords != params_per_line)
- error->all("Incorrect format in COMB potential file");
-
- // words = ptrs to all words in line
-
- nwords = 0;
- words[nwords++] = strtok(line," \t\n\r\f");
- while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
-
- // ielement,jelement,kelement = 1st args
- // if all 3 args are in element list, then parse this line
- // else skip to next line
-
- for (ielement = 0; ielement < nelements; ielement++)
- if (strcmp(words[0],elements[ielement]) == 0) break;
- if (ielement == nelements) continue;
- for (jelement = 0; jelement < nelements; jelement++)
- if (strcmp(words[1],elements[jelement]) == 0) break;
- if (jelement == nelements) continue;
- for (kelement = 0; kelement < nelements; kelement++)
- if (strcmp(words[2],elements[kelement]) == 0) break;
- if (kelement == nelements) continue;
-
- // load up parameter settings and error check their values
-
- if (nparams == maxparam) {
- maxparam += DELTA;
- params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
- "pair:params");
- }
-
- params[nparams].ielement = ielement;
- params[nparams].jelement = jelement;
- params[nparams].kelement = kelement;
- params[nparams].powerm = atof(words[3]);
- params[nparams].c = atof(words[4]);
- params[nparams].d = atof(words[5]);
- params[nparams].h = atof(words[6]);
- params[nparams].powern = atof(words[7]);
- params[nparams].beta = atof(words[8]);
- params[nparams].lam21 = atof(words[9]);
- params[nparams].lam22 = atof(words[10]);
- params[nparams].bigb1 = atof(words[11]);
- params[nparams].bigb2 = atof(words[12]);
- params[nparams].bigr = atof(words[13]);
- params[nparams].bigd = atof(words[14]);
- params[nparams].lam11 = atof(words[15]);
- params[nparams].lam12 = atof(words[16]);
- params[nparams].biga1 = atof(words[17]);
- params[nparams].biga2 = atof(words[18]);
- params[nparams].plp1 = atof(words[19]);
- params[nparams].plp3 = atof(words[20]);
- params[nparams].plp6 = atof(words[21]);
- params[nparams].a123 = atof(words[22]);
- params[nparams].aconf= atof(words[23]);
- params[nparams].addrep = atof(words[24]);
- params[nparams].romigb = atof(words[25]);
- params[nparams].romigc = atof(words[26]);
- params[nparams].romigd = atof(words[27]);
- params[nparams].romiga = atof(words[28]);
- params[nparams].QL1 = atof(words[29]);
- params[nparams].QU1 = atof(words[30]);
- params[nparams].DL1 = atof(words[31]);
- params[nparams].DU1 = atof(words[32]);
- params[nparams].QL2 = atof(words[33]);
- params[nparams].QU2 = atof(words[34]);
- params[nparams].DL2 = atof(words[35]);
- params[nparams].DU2 = atof(words[36]);
- params[nparams].chi = atof(words[37]);
- params[nparams].dj = atof(words[38]);
- params[nparams].dk = atof(words[39]);
- params[nparams].dl = atof(words[40]);
- params[nparams].dm = atof(words[41]);
- params[nparams].esm1 = atof(words[42]);
- params[nparams].cmn1 = atof(words[43]);
- params[nparams].cml1 = atof(words[44]);
- params[nparams].cmn2 = atof(words[45]);
- params[nparams].cml2 = atof(words[46]);
- params[nparams].coulcut = atof(words[47]);
- params[nparams].hfocor = atof(words[48]);
-
- params[nparams].powermint = int(params[nparams].powerm);
-
- // parameter sanity checks
-
- if (params[nparams].lam11 < 0.0 || params[nparams].lam12 < 0.0 ||
- params[nparams].c < 0.0 || params[nparams].d < 0.0 ||
- params[nparams].powern < 0.0 || params[nparams].beta < 0.0 ||
- params[nparams].lam21 < 0.0 || params[nparams].lam22 < 0.0 ||
- params[nparams].bigb1< 0.0 || params[nparams].bigb2< 0.0 ||
- params[nparams].biga1< 0.0 || params[nparams].biga2< 0.0 ||
- params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
- params[nparams].bigd > params[nparams].bigr ||
- params[nparams].powerm - params[nparams].powermint != 0.0 ||
- (params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
- params[nparams].plp1 < 0.0 || params[nparams].plp3 < 0.0 ||
- params[nparams].plp6 < 0.0 ||
- params[nparams].a123 > 360.0 || params[nparams].aconf < 0.0 ||
- params[nparams].addrep < 0.0 || params[nparams].romigb < 0.0 ||
- params[nparams].romigc < 0.0 || params[nparams].romigd < 0.0 ||
- params[nparams].romiga < 0.0 ||
- params[nparams].QL1 > 0.0 || params[nparams].QU1 < 0.0 ||
- params[nparams].DL1 < 0.0 || params[nparams].DU1 > 0.0 ||
- params[nparams].QL2 > 0.0 || params[nparams].QU2 < 0.0 ||
- params[nparams].DL2 < 0.0 || params[nparams].DU2 > 0.0 ||
- params[nparams].chi < 0.0 ||
- params[nparams].dj < 0.0 || params[nparams].dk < 0.0 ||
- params[nparams].dl < 0.0 || params[nparams].dm < 0.0 ||
- params[nparams].esm1 < 0.0)
- error->all("Illegal COMB parameter");
-
- if (params[nparams].lam11 < params[nparams].lam21 ||
- params[nparams].lam12 < params[nparams].lam22 ||
- params[nparams].biga1< params[nparams].bigb1 ||
- params[nparams].biga2< params[nparams].bigb2)
- error->all("Illegal COMB parameter");
-
- nparams++;
- }
-
- delete [] words;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::setup()
-{
- int i,j,k,m,n;
-
- // set elem2param for all element triplet combinations
- // must be a single exact match to lines read from file
- // do not allow for ACB in place of ABC
-
- if (elem2param) memory->destroy_3d_int_array(elem2param);
- elem2param = memory->create_3d_int_array(nelements,nelements,nelements,
- "pair:elem2param");
-
- for (i = 0; i < nelements; i++)
- for (j = 0; j < nelements; j++)
- for (k = 0; k < nelements; k++) {
- n = -1;
- for (m = 0; m < nparams; m++) {
- if (i == params[m].ielement && j == params[m].jelement &&
- k == params[m].kelement) {
- if (n >= 0) error->all("Potential file has duplicate entry");
- n = m;
- }
- }
- if (n < 0) error->all("Potential file is missing an entry");
- elem2param[i][j][k] = n;
- }
-
- // compute parameter values derived from inputs
-
- for (m = 0; m < nparams; m++) {
- params[m].cut = params[m].bigr + params[m].bigd;
- params[m].cutsq = params[m].cut*params[m].cut;
- params[m].c1 = pow(2.0*params[m].powern*1.0e-16,-1.0/params[m].powern);
- params[m].c2 = pow(2.0*params[m].powern*1.0e-8,-1.0/params[m].powern);
- params[m].c3 = 1.0/params[m].c2;
- params[m].c4 = 1.0/params[m].c1;
- params[m].rlm1 = 0.5*(params[m].lam11+params[m].lam12)*params[m].romigc;
- params[m].rlm2 = 0.5*(params[m].lam21+params[m].lam22)*params[m].romigd;
-
- params[m].Qo1 = (params[m].QU1+params[m].QL1)/2.0; // (A22)
- params[m].dQ1 = (params[m].QU1-params[m].QL1)/2.0; // (A21)
- params[m].aB1 = 1.0 /
- (1.0-pow(fabs(params[m].Qo1/params[m].dQ1),10)); // (A20)
- params[m].bB1 = pow(fabs(params[m].aB1),0.1)/params[m].dQ1; // (A19)
- params[m].nD1 = log(params[m].DU1/(params[m].DU1-params[m].DL1))/
- log(params[m].QU1/(params[m].QU1-params[m].QL1));
- params[m].bD1 = (pow((params[m].DL1-params[m].DU1),(1.0/params[m].nD1)))/
- (params[m].QU1-params[m].QL1);
-
- params[m].Qo2 = (params[m].QU2+params[m].QL2)/2.0; // (A22)
- params[m].dQ2 = (params[m].QU2-params[m].QL2)/2.0; // (A21)
- params[m].aB2 = 1.0 /
- (1.0-pow(fabs(params[m].Qo2/params[m].dQ2),10)); // (A20)
- params[m].bB2 = pow(fabs(params[m].aB2),0.1)/params[m].dQ2; // (A19)
- params[m].nD2 = log(params[m].DU2/(params[m].DU2-params[m].DL2))/
- log(params[m].QU2/(params[m].QU2-params[m].QL2));
- params[m].bD2 = (pow((params[m].DL2-params[m].DU2),(1.0/params[m].nD2)))/
- (params[m].QU2-params[m].QL2);
-
- params[m].lcut = params[m].coulcut;
- params[m].lcutsq = params[m].lcut*params[m].lcut;
- }
-
- // set cutmax to max of all params
-
- cutmax = 0.0;
- cor_flag = 0;
- for (m = 0; m < nparams; m++) {
- if (params[m].cut > cutmax) cutmax = params[m].cut;
- if (params[m].lcut > cutmax) cutmax = params[m].lcut;
- if (params[m].hfocor > 0.0001) cor_flag = 1;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::repulsive(Param *param, double rsq, double &fforce,
- int eflag, double &eng, double iq, double jq)
-{
- double r,tmp_fc,tmp_fc_d,tmp_exp,Di,Dj;
- double bigA,Asi,Asj,vrcs,fvrcs,fforce_tmp;
- double rslp,rslp2,rslp4,arr1,arr2,fc2j,fc3j,fcp2j,fcp3j;
-
- double romi = param->addrep;
- double rrcs = param->bigr + param->bigd;
-
- r = sqrt(rsq);
- if (r > rrcs) return ;
-
- tmp_fc = comb_fc(r,param);
- tmp_fc_d = comb_fc_d(r,param);
- tmp_exp = exp(-param->rlm1 * r);
-
- arr1 = 2.22850; arr2 = 1.89350;
- fc2j = comb_fc2(r);
- fc3j = comb_fc3(r);
- fcp2j = comb_fc2_d(r);
- fcp3j = comb_fc3_d(r);
-
- Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-iq)),param->nD1);
- Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-jq)),param->nD2);
- Asi = param->biga1 * exp(param->lam11*Di);
- Asj = param->biga2 * exp(param->lam12*Dj);
-
- if ( Asi > 0.0 && Asj > 0.0 )
- bigA = sqrt(Asi*Asj)*param->romiga;
- else
- bigA = 0.0;
-
- fforce = -bigA * tmp_exp * (tmp_fc_d - tmp_fc*param->rlm1) / r;
-
- // additional repulsion for TiO2 and HfO2 (switch by cor_flag)
-
- vrcs = 0.0; fvrcs = 0.0;
- if (romi > 0.0) {
- if (!cor_flag) {
- vrcs = romi * pow((1.0-r/rrcs),2.0);
- fvrcs= romi * 2.0 * (r/rrcs-1.0)/rrcs; }
- else if (cor_flag) {
- rslp = ((arr1-r)/(arr1-arr2));
- rslp2 = rslp * rslp; rslp4 = rslp2 * rslp2;
- vrcs = fc2j * fc3j * romi * ((50.0*rslp4-30.0*rslp2+4.50))/8.0;
- fvrcs = fcp2j*fcp3j*romi*rslp*(-25.0*rslp2+7.50)/(arr1-arr2);
- }
- fforce_tmp = fforce*vrcs - (tmp_fc * bigA * tmp_exp * fvrcs);
- fforce += fforce_tmp;
- }
-
- // eng = repulsive energy
-
- if (eflag) eng = (tmp_fc * bigA * tmp_exp)*(1.0+vrcs);
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::zeta(Param *param, double rsqij, double rsqik,
- double *delrij, double *delrik)
-{
- double rij,rik,costheta,arg,ex_delr;
-
- rij = sqrt(rsqij);
- if (rij > param->bigr+param->bigd) return 0.0;
- rik = sqrt(rsqik);
- costheta = vec3_dot(delrij,delrik) / (rij*rik);
-
- if (param->powermint == 3) arg = pow(param->rlm2 * (rij-rik),3.0);
- else arg = param->rlm2 * (rij-rik);
-
- if (arg > 69.0776) ex_delr = 1.e30;
- else if (arg < -69.0776) ex_delr = 0.0;
- else ex_delr = exp(arg);
-
- return comb_fc(rik,param) * comb_gijk(costheta,param) * ex_delr;
-}
-
-/* ----------------------------------------------------------------------
- Legendre polynomial bond angle correction to energy
-------------------------------------------------------------------------- */
-
-double PairComb::elp(Param *param, double rsqij, double rsqik,
- double *delrij, double *delrik)
-{
- if (param->aconf > 1.0e-6 || param->plp1 > 1.0e-6 ||
- param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
- double rij,rik,costheta,lp1,lp3,lp6;
- double rmu,rmu2,comtt,fck;
- double pplp1 = param->plp1, pplp3 = param->plp3, pplp6 = param->plp6;
- double c123 = cos(param->a123*PI/180.0);
-
- // cos(theta) of the i-j-k
- // cutoff function of rik
-
- rij = sqrt(rsqij);
- rik = sqrt(rsqik);
- costheta = vec3_dot(delrij,delrik) / (rij*rik);
- fck = comb_fc(rik,param);
- rmu = costheta;
-
- // Legendre Polynomial functions
-
- if (param->plp1 > 1.0e-6 || param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
- rmu2 = rmu*rmu;
- lp1 = rmu; lp3 = 0.5*(5.0*rmu2*rmu-3.0*rmu);
- lp6 = (231.0*rmu2*rmu2*rmu2-315.0*rmu2*rmu2+105.0*rmu2-5.0)/16.0;
- comtt = pplp1*lp1 + pplp3*lp3 + pplp6*lp6;
- } else comtt = 0.0;
-
- // bond-bending terms
-
- if (param->aconf>1e-4) {
- if (param->hfocor >= 0.0)
- comtt += param->aconf *(rmu-c123)*(rmu-c123);
- else if (param->hfocor < 0.0)
- comtt += param->aconf *(4.0-(rmu-c123)*(rmu-c123));
- }
-
- return 0.5 * fck * comtt;
- }
-
- return 0.0;
-}
-
-/* ----------------------------------------------------------------------
- Legendre polynomial bond angle correction to forces
-------------------------------------------------------------------------- */
-
-void PairComb::flp(Param *param, double rsqij, double rsqik,
- double *delrij, double *delrik, double *drilp,
- double *drjlp, double *drklp)
-{
- double ffj1,ffj2,ffk1,ffk2;
- ffj1 = 0.0; ffj2 = 0.0; ffk1 = 0.0; ffk2 = 0.0;
-
- if (param->aconf > 1.0e-4 || param->plp1 > 1.0e-6 ||
- param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
- double rij,rik,costheta,lp1,lp1_d,lp3,lp3_d,lp6,lp6_d;
- double rmu,rmu2,comtt,comtt_d,com4k,com5,fcj,fck,fck_d;
-
- double pplp1 = param->plp1;
- double pplp3 = param->plp3;
- double pplp6 = param->plp6;
- double c123 = cos(param->a123*PI/180.0);
-
- // fck_d = derivative of cutoff function
-
- rij = sqrt(rsqij); rik = sqrt(rsqik);
- costheta = vec3_dot(delrij,delrik) / (rij*rik);
- fcj = comb_fc(rij,param);
- fck = comb_fc(rik,param);
- fck_d = comb_fc_d(rik,param);
- rmu = costheta;
-
- // Legendre Polynomial functions and derivatives
-
- if (param->plp1 > 1.0e-6 || param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
- rmu2 = rmu*rmu;
- lp1 = rmu; lp3 = (2.5*rmu2*rmu-1.5*rmu);
- lp6 = (231.0*rmu2*rmu2*rmu2-315.0*rmu2*rmu2+105.0*rmu2-5.0)/16.0;
- lp1_d = 1.0;lp3_d = (7.5*rmu2-1.5);
- lp6_d = (1386.0*rmu2*rmu2*rmu-1260.0*rmu2*rmu+210.0)/16.0;
- comtt = pplp1*lp1 + pplp3*lp3 + pplp6*lp6;
- comtt_d = pplp1*lp1_d + pplp3*lp3_d + pplp6*lp6_d;
- } else {
- comtt = 0.0;
- comtt_d = 0.0;
- }
-
- // bond-bending terms derivatives
-
- if (param->aconf > 1.0e-4) {
- if (param->hfocor >= 0.0) {
- comtt += param->aconf *(rmu-c123)*(rmu-c123);
- comtt_d += 2.0*param->aconf*(rmu-c123);
- } else if (param->hfocor < 0.0) {
- comtt += param->aconf *(4.0-(rmu-c123)*(rmu-c123));
- comtt_d += -2.0*param->aconf*(rmu-c123);
- }
- }
-
- com4k = fcj * fck_d * comtt;
- com5 = fcj * fck * comtt_d;
-
- ffj1 =-0.5*(com5/(rij*rik));
- ffj2 = 0.5*(com5*rmu/rsqij);
- ffk1 = ffj1;
- ffk2 = 0.5*(-com4k/rik+com5*rmu/rsqik);
-
- } else {
- ffj1 = 0.0; ffj2 = 0.0;
- ffk1 = 0.0; ffk2 = 0.0;
- }
-
- // j-atom
-
- vec3_scale(ffj1,delrik,drjlp); // (k,x[],y[]), y[]=k*x[]
- vec3_scaleadd(ffj2,delrij,drjlp,drjlp); // (k,x[],y[],z[]), z[]=k*x[]+y[]
-
- // k-atom
-
- vec3_scale(ffk1,delrij,drklp);
- vec3_scaleadd(ffk2,delrik,drklp,drklp);
-
- // i-atom
-
- vec3_add(drjlp,drklp,drilp); // (x[],y[],z[]), z[]=x[]+y[]
- vec3_scale(-1.0,drilp,drilp);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::force_zeta(Param *param, double rsq, double zeta_ij,
- double &fforce, double &prefactor,
- int eflag, double &eng, double iq, double jq)
-{
- double r,fa,fa_d,bij;
-
- r = sqrt(rsq);
- if (r > param->bigr+param->bigd) return;
- fa = comb_fa(r,param,iq,jq);
- fa_d = comb_fa_d(r,param,iq,jq);
- bij = comb_bij(zeta_ij,param);
- fforce = 0.5*bij*fa_d / r;
- prefactor = -0.5*fa * comb_bij_d(zeta_ij,param);
-
- // eng = attractive energy
-
- if (eflag) eng = 0.5*bij*fa;
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fc(double r, Param *param)
-{
- double comb_R = param->bigr;
- double comb_D = param->bigd;
-
- if (r < comb_R-comb_D) return 1.0;
- if (r > comb_R+comb_D) return 0.0;
- return 0.5*(1.0 + cos(PI*(r - comb_R)/comb_D));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fc_d(double r, Param *param)
-{
- double comb_R = param->bigr;
- double comb_D = param->bigd;
-
- if (r < comb_R-comb_D) return 0.0;
- if (r > comb_R+comb_D) return 0.0;
- return -(PI2/comb_D) * sin(PI*(r - comb_R)/comb_D);
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fc2(double r)
-{
- double comb_R = 1.89350;
- double comb_D = comb_R + 0.050;
-
- if (r < comb_R) return 0.0;
- if (r > comb_D) return 1.0;
- return 0.5*(1.0 + cos(PI*(r - comb_R)/(comb_D-comb_R)));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fc2_d(double r)
-{
- double comb_R = 1.89350;
- double comb_D = comb_R + 0.050;
-
- if (r < comb_R) return 0.0;
- if (r > comb_D) return 0.0;
- return -(PI2/(comb_D-comb_R)) * sin(PI*(r - comb_R)/(comb_D-comb_R));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fc3(double r)
-{
- double comb_R = 2.51350;
- double comb_D = comb_R + 0.050;
-
- if (r < comb_R) return 1.0;
- if (r > comb_D) return 0.0;
- return 0.5*(1.0 + cos(PI*(r - comb_R)/(comb_D-comb_R)));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fc3_d(double r)
-{
- double comb_R = 2.51350;
- double comb_D = comb_R + 0.050;
-
- if (r < comb_R) return 0.0;
- if (r > comb_D) return 0.0;
- return -(PI2/(comb_D-comb_R)) * sin(PI*(r - comb_R)/(comb_D-comb_R));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::self(Param *param, double qi, double selfpot)
-{
- double self_tmp, cmin, cmax, qmin, qmax;
- double s1=param->chi, s2=param->dj, s3=param->dk, s4=param->dl, s5=param->dm;
-
- self_tmp = 0.0;
- qmin = param->QL1*1.10;
- qmax = param->QU1*0.90;
- cmin = cmax = 3000.0;
-
- self_tmp = qi*(s1+qi*(s2+selfpot+qi*(s3+qi*(s4+qi*qi*s5))));
-
- if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4);
- if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4);
-
- return self_tmp;
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fa(double r, Param *param, double iq, double jq)
-{
- double bigB,Bsi,Bsj;
- double qi,qj,Di,Dj;
-
- if (r > param->bigr + param->bigd) return 0.0;
- qi = iq; qj = jq;
- Di = Dj = Bsi = Bsj = bigB = 0.0;
- Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1);
- Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2);
- Bsi = param->bigb1 * exp(param->lam21*Di)*
- (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10)));
- Bsj = param->bigb2 * exp(param->lam22*Dj)*
- (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10)));
- if (Bsi > 0.0 && Bsj > 0.0) bigB = sqrt(Bsi*Bsj)*param->romigb;
- else bigB = 0.0;
-
- return -bigB * exp(-param->rlm2 * r) * comb_fc(r,param);
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_fa_d(double r, Param *param, double iq, double jq)
-{
- double bigB,Bsi,Bsj;
- double qi,qj,Di,Dj;
-
- if (r > param->bigr + param->bigd) return 0.0;
- qi = iq; qj = jq;
- Di = Dj = Bsi = Bsj = bigB = 0.0;
- Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1);
- Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2);
- Bsi = param->bigb1 * exp(param->lam21*Di)*
- (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10)));
- Bsj = param->bigb2 * exp(param->lam22*Dj)*
- (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10)));
- if (Bsi > 0.0 && Bsj > 0.0) bigB = sqrt(Bsi*Bsj)*param->romigb;
- else bigB = 0.0;
-
- return bigB * exp(-param->rlm2 * r) *
- (param->rlm2 * comb_fc(r,param) - comb_fc_d(r,param));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_bij(double zeta, Param *param)
-{
- double tmp = param->beta * zeta;
- if (tmp > param->c1) return 1.0/sqrt(tmp);
- if (tmp > param->c2)
- return (1.0 - pow(tmp,-param->powern) / (2.0*param->powern))/sqrt(tmp);
- if (tmp < param->c4) return 1.0;
- if (tmp < param->c3)
- return 1.0 - pow(tmp,param->powern)/(2.0*param->powern);
- return pow(1.0 + pow(tmp,param->powern), -1.0/(2.0*param->powern));
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_bij_d(double zeta, Param *param)
-{
- double tmp = param->beta * zeta;
- if (tmp > param->c1) return param->beta * -0.5*pow(tmp,-1.5);
- if (tmp > param->c2)
- return param->beta * (-0.5*pow(tmp,-1.5) *
- (1.0 - 0.5*(1.0 + 1.0/(2.0*param->powern)) *
- pow(tmp,-param->powern)));
- if (tmp < param->c4) return 0.0;
- if (tmp < param->c3)
- return -0.5*param->beta * pow(tmp,param->powern-1.0);
-
- double tmp_n = pow(tmp,param->powern);
- return -0.5 * pow(1.0+tmp_n, -1.0-(1.0/(2.0*param->powern)))*tmp_n / zeta;
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_gijk(double costheta, Param *param)
-{
- double comb_c = param->c;
- double comb_d = param->d;
-
- return (1.0 + pow(comb_c/comb_d,2.0) -
- pow(comb_c,2.0) / (pow(comb_d,2.0) + pow(param->h - costheta,2.0)));
-};
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::comb_gijk_d(double costheta, Param *param)
-{
- double numerator = -2.0 * pow(param->c,2) * (param->h - costheta);
- double denominator = pow(pow(param->d,2.0) +
- pow(param->h - costheta,2.0),2.0);
- return numerator/denominator;
-}
-
-/*------------------------------------------------------------------------- */
-
-void PairComb::attractive(Param *param, double prefactor,
- double rsqij, double rsqik,
- double *delrij, double *delrik,
- double *fi, double *fj, double *fk)
-{
- double rij_hat[3],rik_hat[3];
- double rij,rijinv,rik,rikinv;
-
- rij = sqrt(rsqij);
- rijinv = 1.0/rij;
- vec3_scale(rijinv,delrij,rij_hat);
-
- rik = sqrt(rsqik);
- rikinv = 1.0/rik;
- vec3_scale(rikinv,delrik,rik_hat);
-
- comb_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,param);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::comb_zetaterm_d(double prefactor, double *rij_hat, double rij,
- double *rik_hat, double rik, double *dri,
- double *drj, double *drk, Param *param)
-{
- double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
- double dcosdri[3],dcosdrj[3],dcosdrk[3];
-
- fc = comb_fc(rik,param);
- dfc = comb_fc_d(rik,param);
- if (param->powermint == 3) tmp = pow(param->rlm2 * (rij-rik),3.0);
- else tmp = param->rlm2 * (rij-rik);
-
- if (tmp > 69.0776) ex_delr = 1.e30;
- else if (tmp < -69.0776) ex_delr = 0.0;
- else ex_delr = exp(tmp); // ex_delr is Ygexp
-
- if (param->powermint == 3)
- ex_delr_d = 3.0*pow(param->rlm2,3.0) * pow(rij-rik,2.0)*ex_delr; // com3
- else ex_delr_d = param->rlm2 * ex_delr; // com3
-
- cos_theta = vec3_dot(rij_hat,rik_hat);
- gijk = comb_gijk(cos_theta,param);
- gijk_d = comb_gijk_d(cos_theta,param);
- 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);
- // (k,x[],y[]), y[]=k*x[]
- // (k,x[],y[],z[]), z[]=k*x[]+y[]
-
- 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 PairComb::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);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::sm_table()
-{
- int i,j,k,m,nntypes,ncoul;
- int inty, itype, jtype;
- int iparam_i, iparam_ij, iparam_ji;
- double r,dra,drin,rc,z,zr,zrc,ea,eb,ea3,eb3,alf;
- double exp2er,exp2ersh,fafash,dfafash,F1,dF1,ddF1,E1,E2,E3,E4;
- double exp2ear,exp2ebr,exp2earsh,exp2ebrsh,fafbsh,dfafbsh;
-
- int n = atom->ntypes;
- int *type = atom->type;
-
- dra = 0.001; // lookup table step size
- drin = 0.1; // starting distance of 1/r
- rc = cutmax;
- alf = 0.20;
-
- nntypes = int((n+1)*n/2); // interaction types
- ncoul = int((rc-drin)/dra)+1;
-
- memory->destroy_2d_int_array(intype);
- memory->destroy_2d_double_array(fafb);
- memory->destroy_2d_double_array(dfafb);
- memory->destroy_2d_double_array(ddfafb);
- memory->destroy_2d_double_array(phin);
- memory->destroy_2d_double_array(dphin);
- memory->destroy_2d_double_array(erpaw);
-
- // allocate arrays
-
- intype = memory->create_2d_int_array(n,n,"pair:intype");
- fafb = memory->create_2d_double_array(ncoul,nntypes,"pair:fafb");
- dfafb = memory->create_2d_double_array(ncoul,nntypes,"pair:dfafb");
- ddfafb = memory->create_2d_double_array(ncoul,nntypes,"pair:ddfafb");
- phin = memory->create_2d_double_array(ncoul,nntypes,"pair:phin");
- dphin = memory->create_2d_double_array(ncoul,nntypes,"pair:dphin");
- erpaw = memory->create_2d_double_array(25000,2,"pair:erpaw");
-
- // set interaction number: 0-0=0, 1-1=1, 0-1=1-0=2
-
- m = 0; k = n;
- for (i = 0; i < n; i++) {
- for (j = 0; j < n; j++) {
- if (j == i) {
- intype[i][j] = m;
- m += 1;
- } else if (j != i && j > i) {
- intype[i][j] = k;
- k += 1;
- } else if (j != i && j < i) {
- intype[i][j] = intype[j][i];
- }
- }
- }
-
- // default arrays to zero
-
- for (i = 0; i < ncoul; i ++) {
- for (j = 0; j < nntypes; j ++) {
- fafb[i][j] = 0.0;
- dfafb[i][j] = 0.0;
- ddfafb[i][j] = 0.0;
- phin[i][j] = 0.0;
- dphin[i][j] = 0.0;
- }
- }
-
- // direct 1/r energy with Slater 1S orbital overlap
-
- for (i = 0; i < n; i++) {
- r = drin;
- itype = params[i].ielement;
- iparam_i = elem2param[itype][itype][itype];
- z = params[iparam_i].esm1;
- for (j = 0; j < ncoul; j++) {
- exp2er = exp(-2.0 * z * r);
- phin[j][i] = 1.0 - exp2er * (1.0 + 2.0 * z * r * (1.0 + z * r));
- dphin[j][i] = (4.0 * exp2er * z * z * z * r * r);
- r += dra;
- }
- }
-
- for (i = 0; i < n; i ++) {
- for (j = 0; j < n; j ++) {
- r = drin;
- if (j == i) {
- itype = params[i].ielement;
- inty = intype[itype][itype];
- iparam_i = elem2param[itype][itype][itype];
- z = params[iparam_i].esm1;
- zrc = z * rc;
- exp2ersh = exp(-2.0 * zrc);
- fafash = -exp2ersh * (1.0 / rc +
- z * (11.0/8.0 + 3.0/4.0*zrc + zrc*zrc/6.0));
- dfafash = exp2ersh * (1.0/(rc*rc) + 2.0*z/rc +
- z*z*(2.0 + 7.0/6.0*zrc + zrc*zrc/3.0));
- for (k = 0; k < ncoul; k ++) {
- zr = z * r;
- exp2er = exp(-2.0*zr);
- F1 = -exp2er * (1.0 / r +
- z * (11.0/8.0 + 3.0/4.0*zr + zr*zr/6.0));
- dF1 = exp2er * (1.0/(r*r) + 2.0*z/r +
- z*z*(2.0 + 7.0/6.0*zr + zr*zr/3.0));
- ddF1 = -exp2er * (2.0/(r*r*r) + 4.0*z/(r*r) -
- z*z*z/3.0*(17.0/2.0 + 5.0*zr + 2.0*zr*zr));
- fafb[k][inty] = F1-fafash-(r-rc)*dfafash;
- dfafb[k][inty] = (dF1 - dfafash);
- ddfafb[k][inty] = ddF1;
- r += dra;
- }
- } else if (j != i) {
- itype = params[i].ielement;
- jtype = params[j].ielement;
- inty = intype[itype][jtype];
- iparam_ij = elem2param[itype][jtype][jtype];
- ea = params[iparam_ij].esm1;
- ea3 = ea*ea*ea;
- iparam_ji = elem2param[jtype][itype][itype];
- eb = params[iparam_ji].esm1;
- eb3 = eb*eb*eb;
- E1 = ea*eb3*eb/((ea+eb)*(ea+eb)*(ea-eb)*(ea-eb));
- E2 = eb*ea3*ea/((ea+eb)*(ea+eb)*(eb-ea)*(eb-ea));
- E3 = (3.0*ea*ea*eb3*eb-eb3*eb3) /
- ((ea+eb)*(ea+eb)*(ea+eb)*(ea-eb)*(ea-eb)*(ea-eb));
- E4 = (3.0*eb*eb*ea3*ea-ea3*ea3) /
- ((ea+eb)*(ea+eb)*(ea+eb)*(eb-ea)*(eb-ea)*(eb-ea));
- exp2earsh = exp(-2.0*ea*rc);
- exp2ebrsh = exp(-2.0*eb*rc);
- fafbsh = -exp2earsh*(E1 + E3/rc)-exp2ebrsh*(E2 + E4/rc);
- dfafbsh =
- exp2earsh*(2.0*ea*(E1+E3/rc)+E3/(rc*rc)) +
- exp2ebrsh*(2.0*eb*(E2+E4/rc)+E4/(rc*rc));
- for (k = 0; k < ncoul; k ++) {
- exp2ear = exp(-2.0*ea*r);
- exp2ebr = exp(-2.0*eb*r);
- fafb[k][inty] =
- - exp2ear*(E1+E3/r) - exp2ebr*(E2+E4/r)
- - fafbsh - (r-rc) * dfafbsh;
- dfafb[k][inty] = (exp2ear*(2.0*ea*(E1+E3/r) + E3/(r*r))
- + exp2ebr*(2.0*eb*(E2+E4/r) + E4/(r*r))- dfafbsh);
- ddfafb[k][inty] = (- exp2ear*(E3/(r*r)*(1.0/r+2.0*ea/r+2.0/(r*r))
- + 2.0*ea*(E1+E3/r))-
- exp2ebr*(E4/(r*r)
- *(1.0/r+2.0*eb/r+2.0/(r*r)) +
- 2.0*eb*(E2+E4/r)));
- r += dra;
- }
- }
- }
- }
-
- for (i = 0; i < 25000; i ++) {
- r = dra * i + drin;
- erpaw[i][0] = erfc(r*alf);
- erpaw[i][1] = exp(-r*r*alf*alf);
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::potal_calc(double &calc1, double &calc2, double &calc3)
-{
- double potal,alf,rcoul,fac11,fac11e,esucon;
- int m;
-
- rcoul = 0.0;
- for (m = 0; m < nparams; m++)
- if (params[m].lcut > rcoul) rcoul = params[m].lcut;
-
- alf = 0.20;
- esucon = force->qqr2e;
-
- calc2 = (erfc(rcoul*alf)/rcoul/rcoul+2.0*alf/PIsq*
- exp(-alf*alf*rcoul*rcoul)/rcoul)*esucon/rcoul;
- calc3 = (erfc(rcoul*alf)/rcoul)*esucon;
- calc1 = -(alf/PIsq*esucon+calc3*0.5);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::tri_point(double rsq, int &mr1, int &mr2,
- int &mr3, double &sr1, double &sr2,
- double &sr3, int &itype)
-{
- double r, rin, dr, dd, rr1, rridr, rridr2;
- int m = itype;
-
- rin = 0.10; dr = 0.0010;
- r = sqrt(rsq);
- if (r < rin + 2.0*dr) r = rin + 2.0*dr;
- if (r > cutmax - 2.0*dr) r = cutmax - 2.0*dr;
- rridr = (r-rin)/dr;
-
- mr1 = int(rridr)-1;
- dd = rridr - float(mr1);
- if (dd > 0.5) mr1 += 1;
- mr2 = mr1 + 1;
- mr3 = mr2 + 1;
-
- rr1 = float(mr1)*dr;
- rridr = (r - rin - rr1)/dr;
- rridr2 = rridr * rridr;
-
- sr1 = (rridr2 - rridr) * 0.50;
- sr2 = 1.0 - rridr2;
- sr3 = (rridr2 + rridr) * 0.50;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::direct(int inty, int mr1, int mr2, int mr3, double rsq,
- double sr1, double sr2, double sr3,
- double iq, double jq,
- double potal, double fac11, double fac11e,
- double &pot_tmp, double &pot_d)
-{
- double r,erfcc,fafbn1,potij,sme2,chrij,esucon;
- double r3,erfcd,dfafbn1,smf2,dvdrr,alf,alfdpi;
-
- r = sqrt(rsq);
- r3 = r * rsq;
- alf = 0.20;
- alfdpi = 2.0*alf/PIsq;
- esucon = force->qqr2e;
- pot_tmp = 0.0;
- pot_d = 0.0;
-
- // 1/r energy
-
- erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
- fafbn1= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
- potij = (erfcc/r * esucon - fac11e);
- sme2 = potij + fafbn1 * esucon;
- pot_tmp = 1.0 * iq * jq *sme2;
-
- // 1/r force (wrt r)
-
- erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
- dfafbn1= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
- dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
- smf2 = dvdrr - dfafbn1 * esucon/r;
- pot_d = 1.0 * iq * jq * smf2;
-
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::field(Param *param, double rsq, double iq,double jq,
- double &vionij,double &fvionij)
-{
- double r,r5,r6,rc,rc5,rc6,rf5,drf6,smpn,smpl,rfx1,rfx2;
- double cmi1,cmi2,cmj1,cmj2;
-
- r = sqrt(rsq);
- r5 = r*r*r*r*r;
- r6 = r5 * r;
- rc = param->lcut;
- rc5 = rc*rc*rc*rc*rc;
- rc6 = rc5 * rc;
- cmi1 = param->cmn1;
- cmi2 = param->cmn2;
- cmj1 = param->cml1;
- cmj2 = param->cml2;
- rf5 = 1.0/r5 - 1.0/rc5 + 5.0*(r-rc)/rc6;
- drf6 = 5.0/rc6 - 5.0/r6;
-
- // field correction energy
-
- smpn = rf5*jq*(cmi1+jq*cmi2);
- smpl = rf5*iq*(cmj1+iq*cmj2);
- vionij += 1.0 * (smpn + smpl);
-
- // field correction force
-
- rfx1 = jq*drf6*(cmi1+jq*cmi2)/r;
- rfx2 = iq*drf6*(cmj1+iq*cmj2)/r;
- fvionij -= 1.0 * (rfx1 + rfx2);
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::yasu_char(double *qf_fix, int &igroup)
-{
- int i,j,k,ii,jj,kk,jnum;
- int itag,jtag,itype,jtype,ktype,iparam_i,iparam_ij,iparam_ijk;
- double xtmp,ytmp,ztmp,evdwl,fpair;
- double rsq,rsq1,rsq2,delr1[3],delr2[3],zeta_ij;
- int *ilist,*jlist,*numneigh,**firstneigh;
- double iq,jq,fqi,fqj,fqij,fqjj,ecoul,yaself,yaself_d;
- double potal,fac11,fac11e,sr1,sr2,sr3;
- int mr1,mr2,mr3,inty;
- int zeta_flag;
-
- double **x = atom->x;
- double *q = atom->q;
- int *tag = atom->tag;
- int *type = atom->type;
- int nlocal = atom->nlocal;
-
- int inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- int *mask = atom->mask;
- int groupbit = group->bitmask[igroup];
-
- qf = qf_fix;
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- if (mask[i] & groupbit)
- qf[i] = 0.0;
- }
-
- // communicating charge force to all nodes, first forward then reverse
-
- comm->forward_comm_pair(this);
-
- // self energy correction term: potal
-
- potal_calc(potal,fac11,fac11e);
-
- // loop over full neighbor list of my atoms
-
- fqi = fqj = fqij = fqjj = 0.0;
-
- for (ii = 0; ii < inum; ii ++) {
- i = ilist[ii];
- if (mask[i] & groupbit) {
- itype = map[type[i]];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- iq = q[i];
- iparam_i = elem2param[itype][itype][itype];
-
- // charge force from self energy
-
- fqi = qfo_self(¶ms[iparam_i],iq,potal);
-
- // two-body interactions
-
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
- jtype = map[type[j]];
- jq = q[j];
-
- delr1[0] = x[j][0] - xtmp;
- delr1[1] = x[j][1] - ytmp;
- delr1[2] = x[j][2] - ztmp;
- rsq1 = vec3_dot(delr1,delr1);
-
- iparam_ij = elem2param[itype][jtype][jtype];
-
- // long range q-dependent
-
- if (rsq1 > params[iparam_ij].lcutsq) continue;
-
- inty = intype[itype][jtype];
-
- // polynomial three-point interpolation
-
- tri_point(rsq1,mr1,mr2,mr3,sr1,sr2,sr3,itype);
-
- // 1/r charge forces
-
- qfo_direct(inty,mr1,mr2,mr3,rsq1,sr1,sr2,sr3,fac11e,fqij);
- fqi += jq * fqij; qf[j] += iq * fqij;
-
- // field correction to self energy and charge force
-
- qfo_field(¶ms[iparam_ij],rsq1,iq,jq,fqij,fqjj);
- fqi += fqij;
- qf[j] += fqjj;
-
- // polarization field charge force
- // three-body interactions
-
- if (rsq1 > params[iparam_ij].cutsq) continue;
-
- zeta_ij = 0.0;
-
- for (kk = 0; kk < jnum; kk++) {
- if (jj == kk) continue;
- k = jlist[kk];
- ktype = map[type[k]];
- iparam_ijk = elem2param[itype][jtype][ktype];
-
- delr2[0] = x[k][0] - xtmp;
- delr2[1] = x[k][1] - ytmp;
- delr2[2] = x[k][2] - ztmp;
- rsq2 = vec3_dot(delr2,delr2);
-
- if (rsq2 > params[iparam_ijk].cutsq) continue;
- zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
- }
-
- // charge force in Aij and Bij
-
- qfo_short(¶ms[iparam_ij],rsq1,zeta_ij,iq,jq,fqij,fqjj);
- fqi += fqij; qf[j] += fqjj;
- }
-
- qf[i] += fqi;
-
- }
- }
-
- comm->reverse_comm_pair(this);
-
- // sum charge force on each node and return it
-
- double eneg = 0.0;
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- if (mask[i] & groupbit)
- eneg += qf[i];
- }
- double enegtot;
- MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
- return enegtot;
-}
-
-/* ---------------------------------------------------------------------- */
-
-double PairComb::qfo_self(Param *param, double qi, double selfpot)
-{
- double self_d,cmin,cmax,qmin,qmax;
- double s1 = param->chi;
- double s2 = param->dj;
- double s3 = param->dk;
- double s4 = param->dl;
- double s5 = param->dm;
-
- self_d = 0.0;
- qmin = param->QL1;
- qmax = param->QU1;
- cmin = cmax = 1000.0;
-
- self_d = s1+qi*(2.0*(s2+selfpot)+qi*(3.0*s3+qi*(4.0*s4+qi*qi*6.0*s5)));
-
- /*
- if (qi < qmin) {
- char str[128];
- sprintf(str,"Pair COMB charge %.10f with force %.10f hit min barrier",
- qi,self_d);
- error->warning(str,0);
- self_d += 4.0 * cmin * pow((qi-qmin),3);
- }
- if (qi > qmax) {
- char str[128];
- sprintf(str,"Pair COMB charge %.10f with force %.10f hit max barrier",
- qi,self_d);
- error->warning(str,0);
- self_d += 4.0 * cmax * pow((qi-qmax),3);
- }
- */
-
- return self_d;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::qfo_direct(int inty, int mr1, int mr2, int mr3,
- double rsq, double sr1, double sr2,
- double sr3, double fac11e, double &fqij)
-{
- double r, erfcc, fafbn1, vm, esucon;
-
- r = sqrt(rsq);
- esucon=force->qqr2e;
-
- // 1/r force (wrt q)
-
- erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
- fafbn1= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
- vm = (erfcc/r * esucon - fac11e);
- fqij = 0.5 * (vm+esucon*fafbn1);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::qfo_field(Param *param, double rsq,double iq,double jq,
- double &fqij, double &fqjj)
-{
- double r,r5,r6,rc,rc5,rc6;
- double cmi1,cmi2,cmj1,cmj2,rf5;
-
- fqij = fqjj = 0.0;
- r = sqrt(rsq);
- r5 = r*r*r*r*r;
- r6 = r5 * r;
- rc = param->lcut;
- rc5 = rc*rc*rc*rc*rc;
- rc6 = rc5 * rc;
- cmi1 = param->cmn1;
- cmi2 = param->cmn2;
- cmj1 = param->cml1;
- cmj2 = param->cml2;
- rf5 = 1.0/r5 - 1.0/rc5 + 5.0*(r-rc)/rc6;
-
- // field correction charge force
-
- fqij = 0.5 * rf5 * (cmj1 + 2.0 * iq * cmj2);
- fqjj = 0.5 * rf5 * (cmi1 + 2.0 * jq * cmi2);
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::qfo_short(Param *param, double rsq, double zeta_ij,
- double iq, double jq, double &fqij, double &fqjj)
-{
- double r,tmp_fc,tmp_fc_d,tmp_exp1,tmp_exp2;
- double bigA,Asi,Asj,vrcs;
- double romi = param->addrep,rrcs = param->bigr + param->bigd;
- double qi,qj,Di,Dj,bigB,Bsi,Bsj;
- double QUchi,QOchi,QUchj,QOchj,YYDiqp,YYDjqp;
- double YYAsiqp,YYAsjqp,YYBsiqp,YYBsjqp;
- double caj,cbj,bij,cfqr,cfqs;
- double romie = param->romiga;
- double romib = param->romigb;
- double ca1,ca2,ca3,ca4;
- double rslp,rslp2,rslp4,arr1,arr2,fc2j,fc3j,fcp2j,fcp3j;
-
- qi = iq; qj = jq; r = sqrt(rsq);
- Di = Dj = Asi = Asj = bigA = Bsi = Bsj = bigB = 0.0;
- QUchi = QOchi = QUchj = QOchj = YYDiqp = YYDjqp =0.0;
- YYAsiqp = YYAsjqp = YYBsiqp = YYBsjqp = 0.0;
- caj = cbj = vrcs = cfqr = cfqs = 0.0;
-
- tmp_fc = comb_fc(r,param);
- tmp_fc_d = comb_fc_d(r,param);
- tmp_exp1 = exp(-param->rlm1 * r);
- tmp_exp2 = exp(-param->rlm2 * r);
- bij = comb_bij(zeta_ij,param);
-
- arr1 = 2.22850; arr2 = 1.89350;
- fc2j = comb_fc2(r);
- fc3j = comb_fc3(r);
- fcp2j = comb_fc2_d(r);
- fcp3j = comb_fc3_d(r);
-
- vrcs = 0.0;
- if (romi > 0.0) {
- if (!cor_flag) vrcs = romi * pow((1.0-r/rrcs),2.0);
- else if (cor_flag) {
- rslp = ((arr1-r)/(arr1-arr2));
- rslp2 = rslp * rslp; rslp4 = rslp2 * rslp2;
- vrcs = fc2j * fc3j * romi * ((50.0*rslp4-30.0*rslp2+4.50))/8.0;
- }
- }
-
- Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1);
- Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2);
-
- Asi = param->biga1 * exp(param->lam11*Di);
- Asj = param->biga2 * exp(param->lam12*Dj);
- Bsi = param->bigb1 * exp(param->lam21*Di)*
- (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10)));
- Bsj = param->bigb2 * exp(param->lam22*Dj)*
- (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10)));
-
- QUchi = (param->QU1-qi)*param->bD1;
- QUchj = (param->QU2-qj)*param->bD2;
- QOchi = (qi-param->Qo1)*param->bB1;
- QOchj = (qj-param->Qo2)*param->bB2;
-
- if (QUchi == 0.0) YYDiqp = 0.0;
- else YYDiqp = -param->nD1 * QUchi * param->bD1 *
- pow(fabs(QUchi),(param->nD1-2.0));
-
- if (QUchj == 0.0) YYDjqp = 0.0;
- else YYDjqp = -param->nD2 * QUchj * param->bD2 *
- pow(fabs(QUchj),(param->nD2-2.0));
-
- YYAsiqp = Asi * param->lam11 * YYDiqp;
- YYAsjqp = Asj * param->lam12 * YYDjqp;
-
- if (QOchi == 0.0)
- YYBsiqp=Bsi*param->lam21*YYDiqp;
- else
- YYBsiqp=Bsi*param->lam21*YYDiqp-param->bigb1*exp(param->lam21*Di)*
- 10.0*QOchi*param->bB1*pow(fabs(QOchi),(10.0-2.0));
-
- if (QOchj == 0.0)
- YYBsjqp=Bsj*param->lam22*YYDjqp;
- else
- YYBsjqp=Bsj*param->lam22*YYDjqp-param->bigb2*exp(param->lam22*Dj)*
- 10.0*QOchj*param->bB2*pow(fabs(QOchj),(10.0-2.0));
-
- if (Asi > 0.0 && Asj > 0.0) caj = 1.0/(2.0*sqrt(Asi*Asj)) * romie;
- else caj == 0.0;
-
- if (Bsi > 0.0 && Bsj > 0.0) cbj = 1.0/(2.0*sqrt(Bsi*Bsj)) * romib ;
- else cbj == 0.0;
-
- cfqr = 0.50 * tmp_fc * (1.0 + vrcs); // 0.5 b/c full atom loop
- cfqs = -0.50 * tmp_fc * bij;
-
- ca1 = Asj * caj * YYAsiqp;
- ca2 = Bsj * cbj * YYBsiqp;
- ca3 = Asi * caj * YYAsjqp;
- ca4 = Bsi * cbj * YYBsjqp;
-
- fqij = cfqr * tmp_exp1 * ca1;
- fqij += cfqs * tmp_exp2 * ca2;
- fqjj = cfqr * tmp_exp1 * ca3;
- fqjj += cfqs * tmp_exp2 * ca4;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::Over_cor(Param *param, double rsq1, int NCoi,
- double &Eov, double &Fov)
-{
- double ECo,BCo,tmp_fc,tmp_fc_d;
- double r = sqrt(rsq1);
- int NCon = NCoi - 7;
-
- tmp_fc = comb_fc(r,param);
- tmp_fc_d = comb_fc(r,param);
- Eov = 0.0; Fov = 0.0;
- ECo = param->hfocor;
- BCo = 0.1;
-
- if (NCon >= 0.20) {
- Eov = tmp_fc * ECo * NCon/(1.0+exp(BCo*NCon));
- Fov = -(tmp_fc_d*Eov + tmp_fc*ECo/(1.0+exp(BCo*NCon)) -
- (tmp_fc*ECo*NCon*BCo*exp(BCo*NCon)) /
- ((1.0+exp(BCo*NCon))*(1.0+exp(BCo*NCon))));
- Fov /= r;
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int PairComb::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
-{
- int i,j,m;
-
- m = 0;
- for (i = 0; i < n; i ++) {
- j = list[i];
- buf[m++] = qf[j];
- }
- return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::unpack_comm(int n, int first, double *buf)
-{
- int i,m,last;
-
- m = 0;
- last = first + n ;
- for (i = first; i < last; i++) qf[i] = buf[m++];
-}
-
-/* ---------------------------------------------------------------------- */
-
-int PairComb::pack_reverse_comm(int n, int first, double *buf)
-{
- int i,m,last;
-
- m = 0;
- last = first + n;
- for (i = first; i < last; i++) buf[m++] = qf[i];
- return 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairComb::unpack_reverse_comm(int n, int *list, double *buf)
-{
- int i,j,m;
-
- m = 0;
- for (i = 0; i < n; i++) {
- j = list[i];
- qf[j] += buf[m++];
- }
-}
-
-/* ----------------------------------------------------------------------
- memory usage of local atom-based arrays
-------------------------------------------------------------------------- */
-
-double PairComb::memory_usage()
-{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
- bytes += nmax * sizeof(int);
- return bytes;
-}
+/* ----------------------------------------------------------------------
+ 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: Tzu-Ray Shan (U Florida, rayshan@ufl.edu)
+ LAMMPS implementation of the Charge-optimized many-body (COMB) potential
+ based on the HELL MD program (Prof Simon Phillpot, UF, sphil@mse.ufl.edu)
+ and Aidan Thompson's Tersoff code in LAMMPS
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_comb.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "memory.h"
+#include "error.h"
+#include "group.h"
+
+using namespace LAMMPS_NS;
+
+#define MAXLINE 1024
+#define DELTA 4
+
+/* ---------------------------------------------------------------------- */
+
+PairComb::PairComb(LAMMPS *lmp) : Pair(lmp)
+{
+ single_enable = 0;
+ one_coeff = 1;
+
+ PI = 4.0*atan(1.0);
+ PI2 = 2.0*atan(1.0);
+ PI4 = atan(1.0);
+ PIsq = sqrt(PI);
+
+ nmax = 0;
+ NCo = NULL;
+
+ nelements = 0;
+ elements = NULL;
+ nparams = 0;
+ maxparam = 0;
+ params = NULL;
+ elem2param = NULL;
+
+ intype = NULL;
+ fafb = NULL;
+ dfafb = NULL;
+ ddfafb = NULL;
+ phin = NULL;
+ dphin = NULL;
+ erpaw = NULL;
+
+ // set comm size needed by this Pair
+
+ comm_forward = 1;
+ comm_reverse = 1;
+}
+
+/* ----------------------------------------------------------------------
+ check if allocated, since class can be destructed when incomplete
+------------------------------------------------------------------------- */
+
+PairComb::~PairComb()
+{
+ memory->sfree(NCo);
+
+ if (elements)
+ for (int i = 0; i < nelements; i++) delete [] elements[i];
+ delete [] elements;
+ memory->sfree(params);
+ memory->destroy_3d_int_array(elem2param);
+
+ memory->destroy_2d_int_array(intype);
+ memory->destroy_2d_double_array(fafb);
+ memory->destroy_2d_double_array(dfafb);
+ memory->destroy_2d_double_array(ddfafb);
+ memory->destroy_2d_double_array(phin);
+ memory->destroy_2d_double_array(dphin);
+ memory->destroy_2d_double_array(erpaw);
+
+ if (allocated) {
+ memory->destroy_2d_int_array(setflag);
+ memory->destroy_2d_double_array(cutsq);
+ delete [] map;
+ delete [] esm;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::compute(int eflag, int vflag)
+{
+ int i,j,k,ii,jj,kk,inum,jnum,iparam_i;
+ int itag,jtag,itype,jtype,ktype,iparam_ij,iparam_ijk;
+ double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair;
+ double rsq,rsq1,rsq2;
+ double delr1[3],delr2[3],fi[3],fj[3],fk[3];
+ double zeta_ij,prefactor;
+ double fqi,fqij;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ int mr1,mr2,mr3;
+ int rsc,inty;
+ double elp_ij,filp[3],fjlp[3],fklp[3];
+ double iq,jq;
+ double yaself;
+ double potal,fac11,fac11e;
+ double vionij,fvionij,sr1,sr2,sr3,Eov,Fov;
+
+ evdwl = ecoul = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = vflag_fdotr = vflag_atom = 0;
+
+ // grow coordination array if necessary
+
+ if (atom->nmax > nmax) {
+ memory->sfree(NCo);
+ nmax = atom->nmax;
+ NCo = (int *) memory->smalloc(nmax*sizeof(double),"pair:NCo");
+ }
+
+ double **x = atom->x;
+ double **f = atom->f;
+ double *q = atom->q;
+ int *tag = atom->tag;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+ int newton_pair = force->newton_pair;
+ int ntype = atom->ntypes;
+
+ inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ yaself = vionij = fvionij = Eov = Fov = 0.0;
+
+ // self energy correction term: potal
+
+ potal_calc(potal,fac11,fac11e);
+
+ // 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];
+ iq = q[i];
+ NCo[i] = 0;
+ iparam_i = elem2param[itype][itype][itype];
+
+ // self energy, only on i atom
+
+ yaself = self(¶ms[iparam_i],iq,potal);
+
+ if (evflag) ev_tally(i,i,nlocal,0,yaself,0.0,0.0,0.0,0.0,0.0);
+
+ // two-body interactions (long and short repulsive)
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ 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;
+ }
+
+ // Qj calculates 2-body Coulombic
+
+ jtype = map[type[j]];
+ jq = q[j];
+
+ 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[itype][jtype][jtype];
+
+ // long range q-dependent
+
+ if (rsq > params[iparam_ij].lcutsq) continue;
+
+ inty = intype[itype][jtype];
+
+ // polynomial three-point interpolation
+
+ tri_point(rsq, mr1, mr2, mr3, sr1, sr2, sr3, itype);
+
+ // 1/r energy and forces
+
+ direct(inty,mr1,mr2,mr3,rsq,sr1,sr2,sr3,iq,jq,
+ potal,fac11,fac11e,vionij,fvionij);
+
+ // field correction to self energy
+
+ field(¶ms[iparam_ij],rsq,iq,jq,vionij,fvionij);
+
+ // polarization field
+ // sums up long range forces
+
+ f[i][0] += delx*fvionij;
+ f[i][1] += dely*fvionij;
+ f[i][2] += delz*fvionij;
+ f[j][0] -= delx*fvionij;
+ f[j][1] -= dely*fvionij;
+ f[j][2] -= delz*fvionij;
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,0.0,vionij,fvionij,delx,dely,delz);
+
+ // short range q-independent
+
+ if (rsq > params[iparam_ij].cutsq) continue;
+
+ repulsive(¶ms[iparam_ij],rsq,fpair,eflag,evdwl,iq,jq);
+
+ // repulsion is pure two-body, sums up pair repulsive forces
+
+ 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);
+ }
+
+ // accumulate coordination number information
+
+ if (cor_flag) {
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ jtype = map[type[j]];
+ iparam_ij = elem2param[itype][jtype][jtype];
+
+ if(params[iparam_ij].hfocor > 0.0 ) {
+ delr1[0] = x[j][0] - xtmp;
+ delr1[1] = x[j][1] - ytmp;
+ delr1[2] = x[j][2] - ztmp;
+ rsq1 = vec3_dot(delr1,delr1);
+
+ if (rsq1 > params[iparam_ij].cutsq) continue;
+ NCo[i] += 1;
+ }
+ }
+ }
+
+ // three-body interactions
+ // skip immediately if I-J is not within cutoff
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ jtype = map[type[j]];
+ iparam_ij = elem2param[itype][jtype][jtype];
+
+ // this Qj for q-dependent BSi
+
+ jq = q[j];
+
+ delr1[0] = x[j][0] - xtmp;
+ delr1[1] = x[j][1] - ytmp;
+ delr1[2] = x[j][2] - ztmp;
+ rsq1 = vec3_dot(delr1,delr1);
+
+ if (rsq1 > params[iparam_ij].cutsq) continue;
+
+ // accumulate bondorder zeta for each i-j interaction via loop over k
+
+ zeta_ij = 0.0;
+ cuo_flag1 = 0; cuo_flag2 = 0;
+
+ for (kk = 0; kk < jnum; kk++) {
+ if (jj == kk) continue;
+ k = jlist[kk];
+ ktype = map[type[k]];
+ iparam_ijk = elem2param[itype][jtype][ktype];
+
+ delr2[0] = x[k][0] - xtmp;
+ delr2[1] = x[k][1] - ytmp;
+ delr2[2] = x[k][2] - ztmp;
+ rsq2 = vec3_dot(delr2,delr2);
+
+ if (rsq2 > params[iparam_ijk].cutsq) continue;
+
+ zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
+
+ if (params[iparam_ijk].hfocor == -2.0) cuo_flag1 = 1;
+ if (params[iparam_ijk].hfocor == -1.0) cuo_flag2 = 1;
+ }
+
+ if (cuo_flag1 && cuo_flag2) cuo_flag = 1;
+ else cuo_flag = 0;
+
+ force_zeta(¶ms[iparam_ij],rsq1,zeta_ij,fpair,
+ prefactor,eflag,evdwl,iq,jq);
+
+ // over-coordination correction for HfO2
+
+ if (cor_flag && NCo[i] != 0)
+ Over_cor(¶ms[iparam_ij],rsq1,NCo[i],Eov, Fov);
+ evdwl += Eov;
+ fpair += Fov;
+
+ 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 (3-body forces)
+
+ for (kk = 0; kk < jnum; kk++) {
+ if (jj == kk) continue;
+ k = jlist[kk];
+ ktype = map[type[k]];
+ iparam_ijk = elem2param[itype][jtype][ktype];
+
+ delr2[0] = x[k][0] - xtmp;
+ delr2[1] = x[k][1] - ytmp;
+ delr2[2] = x[k][2] - ztmp;
+ rsq2 = vec3_dot(delr2,delr2);
+ if (rsq2 > params[iparam_ijk].cutsq) continue;
+
+ for (rsc = 0; rsc < 3; rsc++)
+ fi[rsc] = fj[rsc] = fk[rsc] = 0.0;
+
+ attractive(¶ms[iparam_ijk],prefactor,
+ rsq1,rsq2,delr1,delr2,fi,fj,fk);
+
+ // 3-body LP and BB correction and forces
+
+ elp_ij = elp(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
+ flp(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2,filp,fjlp,fklp);
+
+ for (rsc = 0; rsc < 3; rsc++) {
+ fi[rsc] += filp[rsc];
+ fj[rsc] += fjlp[rsc];
+ fk[rsc] += fklp[rsc];
+ }
+
+ for (rsc = 0; rsc < 3; rsc++) {
+ f[i][rsc] += fi[rsc];
+ f[j][rsc] += fj[rsc];
+ f[k][rsc] += fk[rsc];
+ }
+
+ if (evflag)
+ ev_tally(i,j,nlocal,newton_pair,elp_ij,0.0,0.0,0.0,0.0,0.0);
+ if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);
+
+ }
+ }
+
+ if (cuo_flag) params[iparam_i].cutsq *= 0.65;
+ }
+
+ cuo_flag = 0;
+
+ if (vflag_fdotr) virial_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::allocate()
+{
+ allocated = 1;
+ int n = atom->ntypes;
+
+ setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
+ cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+
+ map = new int[n+1];
+ esm = new double[n];
+}
+
+/* ----------------------------------------------------------------------
+ global settings
+------------------------------------------------------------------------- */
+
+void PairComb::settings(int narg, char **arg)
+{
+ if (narg > 0) error->all("Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairComb::coeff(int narg, char **arg)
+{
+ int i,j,n;
+ double r;
+
+ if (!allocated) allocate();
+
+ if (narg != 3 + atom->ntypes)
+ error->all("Incorrect args for pair coefficients");
+
+ // insure I,J args are * *
+
+ if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
+ error->all("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();
+
+ n = atom->ntypes;
+
+ // generate streitz-mintmire direct 1/r energy look-up table
+
+ if (comm->me == 0 && screen) fprintf(screen,"Pair COMB:\n");
+ if (comm->me == 0 && screen)
+ fprintf(screen," generating Coulomb integral lookup table ...\n");
+ sm_table();
+
+ if (cor_flag && comm->me == 0 && screen)
+ fprintf(screen," will apply over-coordination correction ...\n");
+ if (!cor_flag&& comm->me == 0 && screen)
+ fprintf(screen," will not apply over-coordination correction ...\n");
+
+ // clear setflag since coeff() called once with I,J = * *
+
+ 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("Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+ init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairComb::init_style()
+{
+ if (atom->tag_enable == 0)
+ error->all("Pair style COMB requires atom IDs");
+ if (force->newton_pair == 0)
+ error->all("Pair style COMB requires newton pair on");
+ if (!atom->q_flag)
+ error->all("Pair style COMB requires atom attribute q");
+
+ // ptr to QEQ fix
+
+ //for (i = 0; i < modify->nfix; i++)
+ // if (strcmp(modify->fix[i]->style,"qeq") == 0) break;
+ //if (i < modify->nfix) fixqeq = (FixQEQ *) modify->fix[i];
+ //else fixqeq = NULL;
+
+ // 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 PairComb::init_one(int i, int j)
+{
+ if (setflag[i][j] == 0) error->all("All pair coeffs are not set");
+ return cutmax;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::read_file(char *file)
+{
+ int params_per_line = 49;
+ char **words = new char*[params_per_line+1];
+
+ if (params) delete [] params;
+ params = NULL;
+ nparams = 0;
+
+ // open file on proc 0
+
+ FILE *fp;
+ if (comm->me == 0) {
+ fp = fopen(file,"r");
+ if (fp == NULL) {
+ char str[128];
+ sprintf(str,"Cannot open COMB potential file %s",file);
+ error->one(str);
+ }
+ }
+
+ // read each line out of file, skipping blank lines or leading '#'
+ // store line of params if all 3 element tags are in element list
+
+ int n,nwords,ielement,jelement,kelement;
+ char line[MAXLINE],*ptr;
+ int eof = 0;
+
+ while (1) {
+ if (comm->me == 0) {
+ ptr = fgets(line,MAXLINE,fp);
+ if (ptr == NULL) {
+ eof = 1;
+ fclose(fp);
+ } else n = strlen(line) + 1;
+ }
+ MPI_Bcast(&eof,1,MPI_INT,0,world);
+ if (eof) break;
+ MPI_Bcast(&n,1,MPI_INT,0,world);
+ MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+ // strip comment, skip line if blank
+
+ if (ptr = strchr(line,'#')) *ptr = '\0';
+ nwords = atom->count_words(line);
+ if (nwords == 0) continue;
+
+ // concatenate additional lines until have params_per_line words
+
+ while (nwords < params_per_line) {
+ n = strlen(line);
+ if (comm->me == 0) {
+ ptr = fgets(&line[n],MAXLINE-n,fp);
+ if (ptr == NULL) {
+ eof = 1;
+ fclose(fp);
+ } else n = strlen(line) + 1;
+ }
+ MPI_Bcast(&eof,1,MPI_INT,0,world);
+ if (eof) break;
+ MPI_Bcast(&n,1,MPI_INT,0,world);
+ MPI_Bcast(line,n,MPI_CHAR,0,world);
+ if (ptr = strchr(line,'#')) *ptr = '\0';
+ nwords = atom->count_words(line);
+ }
+
+ if (nwords != params_per_line)
+ error->all("Incorrect format in COMB potential file");
+
+ // words = ptrs to all words in line
+
+ nwords = 0;
+ words[nwords++] = strtok(line," \t\n\r\f");
+ while (words[nwords++] = strtok(NULL," \t\n\r\f")) continue;
+
+ // ielement,jelement,kelement = 1st args
+ // if all 3 args are in element list, then parse this line
+ // else skip to next line
+
+ for (ielement = 0; ielement < nelements; ielement++)
+ if (strcmp(words[0],elements[ielement]) == 0) break;
+ if (ielement == nelements) continue;
+ for (jelement = 0; jelement < nelements; jelement++)
+ if (strcmp(words[1],elements[jelement]) == 0) break;
+ if (jelement == nelements) continue;
+ for (kelement = 0; kelement < nelements; kelement++)
+ if (strcmp(words[2],elements[kelement]) == 0) break;
+ if (kelement == nelements) continue;
+
+ // load up parameter settings and error check their values
+
+ if (nparams == maxparam) {
+ maxparam += DELTA;
+ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+ "pair:params");
+ }
+
+ params[nparams].ielement = ielement;
+ params[nparams].jelement = jelement;
+ params[nparams].kelement = kelement;
+ params[nparams].powerm = atof(words[3]);
+ params[nparams].c = atof(words[4]);
+ params[nparams].d = atof(words[5]);
+ params[nparams].h = atof(words[6]);
+ params[nparams].powern = atof(words[7]);
+ params[nparams].beta = atof(words[8]);
+ params[nparams].lam21 = atof(words[9]);
+ params[nparams].lam22 = atof(words[10]);
+ params[nparams].bigb1 = atof(words[11]);
+ params[nparams].bigb2 = atof(words[12]);
+ params[nparams].bigr = atof(words[13]);
+ params[nparams].bigd = atof(words[14]);
+ params[nparams].lam11 = atof(words[15]);
+ params[nparams].lam12 = atof(words[16]);
+ params[nparams].biga1 = atof(words[17]);
+ params[nparams].biga2 = atof(words[18]);
+ params[nparams].plp1 = atof(words[19]);
+ params[nparams].plp3 = atof(words[20]);
+ params[nparams].plp6 = atof(words[21]);
+ params[nparams].a123 = atof(words[22]);
+ params[nparams].aconf= atof(words[23]);
+ params[nparams].addrep = atof(words[24]);
+ params[nparams].romigb = atof(words[25]);
+ params[nparams].romigc = atof(words[26]);
+ params[nparams].romigd = atof(words[27]);
+ params[nparams].romiga = atof(words[28]);
+ params[nparams].QL1 = atof(words[29]);
+ params[nparams].QU1 = atof(words[30]);
+ params[nparams].DL1 = atof(words[31]);
+ params[nparams].DU1 = atof(words[32]);
+ params[nparams].QL2 = atof(words[33]);
+ params[nparams].QU2 = atof(words[34]);
+ params[nparams].DL2 = atof(words[35]);
+ params[nparams].DU2 = atof(words[36]);
+ params[nparams].chi = atof(words[37]);
+ params[nparams].dj = atof(words[38]);
+ params[nparams].dk = atof(words[39]);
+ params[nparams].dl = atof(words[40]);
+ params[nparams].dm = atof(words[41]);
+ params[nparams].esm1 = atof(words[42]);
+ params[nparams].cmn1 = atof(words[43]);
+ params[nparams].cml1 = atof(words[44]);
+ params[nparams].cmn2 = atof(words[45]);
+ params[nparams].cml2 = atof(words[46]);
+ params[nparams].coulcut = atof(words[47]);
+ params[nparams].hfocor = atof(words[48]);
+
+ params[nparams].powermint = int(params[nparams].powerm);
+
+ // parameter sanity checks
+
+ if (params[nparams].lam11 < 0.0 || params[nparams].lam12 < 0.0 ||
+ params[nparams].c < 0.0 || params[nparams].d < 0.0 ||
+ params[nparams].powern < 0.0 || params[nparams].beta < 0.0 ||
+ params[nparams].lam21 < 0.0 || params[nparams].lam22 < 0.0 ||
+ params[nparams].bigb1< 0.0 || params[nparams].bigb2< 0.0 ||
+ params[nparams].biga1< 0.0 || params[nparams].biga2< 0.0 ||
+ params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
+ params[nparams].bigd > params[nparams].bigr ||
+ params[nparams].powerm - params[nparams].powermint != 0.0 ||
+ (params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
+ params[nparams].plp1 < 0.0 || params[nparams].plp3 < 0.0 ||
+ params[nparams].plp6 < 0.0 ||
+ params[nparams].a123 > 360.0 || params[nparams].aconf < 0.0 ||
+ params[nparams].addrep < 0.0 || params[nparams].romigb < 0.0 ||
+ params[nparams].romigc < 0.0 || params[nparams].romigd < 0.0 ||
+ params[nparams].romiga < 0.0 ||
+ params[nparams].QL1 > 0.0 || params[nparams].QU1 < 0.0 ||
+ params[nparams].DL1 < 0.0 || params[nparams].DU1 > 0.0 ||
+ params[nparams].QL2 > 0.0 || params[nparams].QU2 < 0.0 ||
+ params[nparams].DL2 < 0.0 || params[nparams].DU2 > 0.0 ||
+ params[nparams].chi < 0.0 ||
+// params[nparams].dj < 0.0 || params[nparams].dk < 0.0 ||
+// params[nparams].dl < 0.0 || params[nparams].dm < 0.0 ||
+ params[nparams].esm1 < 0.0)
+ error->all("Illegal COMB parameter");
+
+ if (params[nparams].lam11 < params[nparams].lam21 ||
+ params[nparams].lam12 < params[nparams].lam22 ||
+ params[nparams].biga1< params[nparams].bigb1 ||
+ params[nparams].biga2< params[nparams].bigb2)
+ error->all("Illegal COMB parameter");
+
+ nparams++;
+ }
+
+ delete [] words;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::setup()
+{
+ int i,j,k,m,n;
+
+ // set elem2param for all element triplet combinations
+ // must be a single exact match to lines read from file
+ // do not allow for ACB in place of ABC
+
+ if (elem2param) memory->destroy_3d_int_array(elem2param);
+ elem2param = memory->create_3d_int_array(nelements,nelements,nelements,
+ "pair:elem2param");
+
+ for (i = 0; i < nelements; i++)
+ for (j = 0; j < nelements; j++)
+ for (k = 0; k < nelements; k++) {
+ n = -1;
+ for (m = 0; m < nparams; m++) {
+ if (i == params[m].ielement && j == params[m].jelement &&
+ k == params[m].kelement) {
+ if (n >= 0) error->all("Potential file has duplicate entry");
+ n = m;
+ }
+ }
+ if (n < 0) error->all("Potential file is missing an entry");
+ elem2param[i][j][k] = n;
+ }
+
+ // compute parameter values derived from inputs
+
+ for (m = 0; m < nparams; m++) {
+ params[m].cut = params[m].bigr + params[m].bigd;
+ params[m].cutsq = params[m].cut*params[m].cut;
+ params[m].c1 = pow(2.0*params[m].powern*1.0e-16,-1.0/params[m].powern);
+ params[m].c2 = pow(2.0*params[m].powern*1.0e-8,-1.0/params[m].powern);
+ params[m].c3 = 1.0/params[m].c2;
+ params[m].c4 = 1.0/params[m].c1;
+ params[m].rlm1 = 0.5*(params[m].lam11+params[m].lam12)*params[m].romigc;
+ params[m].rlm2 = 0.5*(params[m].lam21+params[m].lam22)*params[m].romigd;
+
+ params[m].Qo1 = (params[m].QU1+params[m].QL1)/2.0; // (A22)
+ params[m].dQ1 = (params[m].QU1-params[m].QL1)/2.0; // (A21)
+ params[m].aB1 = 1.0 /
+ (1.0-pow(fabs(params[m].Qo1/params[m].dQ1),10)); // (A20)
+ params[m].bB1 = pow(fabs(params[m].aB1),0.1)/params[m].dQ1; // (A19)
+ params[m].nD1 = log(params[m].DU1/(params[m].DU1-params[m].DL1))/
+ log(params[m].QU1/(params[m].QU1-params[m].QL1));
+ params[m].bD1 = (pow((params[m].DL1-params[m].DU1),(1.0/params[m].nD1)))/
+ (params[m].QU1-params[m].QL1);
+
+ params[m].Qo2 = (params[m].QU2+params[m].QL2)/2.0; // (A22)
+ params[m].dQ2 = (params[m].QU2-params[m].QL2)/2.0; // (A21)
+ params[m].aB2 = 1.0 /
+ (1.0-pow(fabs(params[m].Qo2/params[m].dQ2),10)); // (A20)
+ params[m].bB2 = pow(fabs(params[m].aB2),0.1)/params[m].dQ2; // (A19)
+ params[m].nD2 = log(params[m].DU2/(params[m].DU2-params[m].DL2))/
+ log(params[m].QU2/(params[m].QU2-params[m].QL2));
+ params[m].bD2 = (pow((params[m].DL2-params[m].DU2),(1.0/params[m].nD2)))/
+ (params[m].QU2-params[m].QL2);
+
+ params[m].lcut = params[m].coulcut;
+ params[m].lcutsq = params[m].lcut*params[m].lcut;
+ }
+
+ // set cutmax to max of all params
+
+ cutmax = 0.0;
+ cor_flag = 0;
+ for (m = 0; m < nparams; m++) {
+ if (params[m].cut > cutmax) cutmax = params[m].cut;
+ if (params[m].lcut > cutmax) cutmax = params[m].lcut;
+ if (params[m].hfocor > 0.0001) cor_flag = 1;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::repulsive(Param *param, double rsq, double &fforce,
+ int eflag, double &eng, double iq, double jq)
+{
+ double r,tmp_fc,tmp_fc_d,tmp_exp,Di,Dj;
+ double bigA,Asi,Asj,vrcs,fvrcs,fforce_tmp;
+ double rslp,rslp2,rslp4,arr1,arr2,fc2j,fc3j,fcp2j,fcp3j;
+
+ double romi = param->addrep;
+ double rrcs = param->bigr + param->bigd;
+
+ r = sqrt(rsq);
+ if (r > rrcs) return ;
+
+ tmp_fc = comb_fc(r,param);
+ tmp_fc_d = comb_fc_d(r,param);
+ tmp_exp = exp(-param->rlm1 * r);
+
+ arr1 = 2.22850; arr2 = 1.89350;
+ fc2j = comb_fc2(r);
+ fc3j = comb_fc3(r);
+ fcp2j = comb_fc2_d(r);
+ fcp3j = comb_fc3_d(r);
+
+ Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-iq)),param->nD1);
+ Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-jq)),param->nD2);
+ Asi = param->biga1 * exp(param->lam11*Di);
+ Asj = param->biga2 * exp(param->lam12*Dj);
+
+ if ( Asi > 0.0 && Asj > 0.0 )
+ bigA = sqrt(Asi*Asj)*param->romiga;
+ else
+ bigA = 0.0;
+
+ fforce = -bigA * tmp_exp * (tmp_fc_d - tmp_fc*param->rlm1) / r;
+
+ // additional repulsion for TiO2 and HfO2 (switch by cor_flag)
+
+ vrcs = 0.0; fvrcs = 0.0;
+ if (romi > 0.0) {
+ if (!cor_flag) {
+ vrcs = romi * pow((1.0-r/rrcs),2.0);
+ fvrcs= romi * 2.0 * (r/rrcs-1.0)/rrcs; }
+ else if (cor_flag) {
+ rslp = ((arr1-r)/(arr1-arr2));
+ rslp2 = rslp * rslp; rslp4 = rslp2 * rslp2;
+ vrcs = fc2j * fc3j * romi * ((50.0*rslp4-30.0*rslp2+4.50))/8.0;
+ fvrcs = fcp2j*fcp3j*romi*rslp*(-25.0*rslp2+7.50)/(arr1-arr2);
+ }
+ fforce_tmp = fforce*vrcs - (tmp_fc * bigA * tmp_exp * fvrcs);
+ fforce += fforce_tmp;
+ }
+
+ // eng = repulsive energy
+
+ if (eflag) eng = (tmp_fc * bigA * tmp_exp)*(1.0+vrcs);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::zeta(Param *param, double rsqij, double rsqik,
+ double *delrij, double *delrik)
+{
+ double rij,rik,costheta,arg,ex_delr;
+
+ rij = sqrt(rsqij);
+ if (rij > param->bigr+param->bigd) return 0.0;
+ rik = sqrt(rsqik);
+ costheta = vec3_dot(delrij,delrik) / (rij*rik);
+
+ if (param->powermint == 3) arg = pow(param->rlm2 * (rij-rik),3.0);
+ else arg = param->rlm2 * (rij-rik);
+
+ if (arg > 69.0776) ex_delr = 1.e30;
+ else if (arg < -69.0776) ex_delr = 0.0;
+ else ex_delr = exp(arg);
+
+ return comb_fc(rik,param) * comb_gijk(costheta,param) * ex_delr;
+}
+
+/* ----------------------------------------------------------------------
+ Legendre polynomial bond angle correction to energy
+------------------------------------------------------------------------- */
+
+double PairComb::elp(Param *param, double rsqij, double rsqik,
+ double *delrij, double *delrik)
+{
+ if (param->aconf > 1.0e-6 || param->plp1 > 1.0e-6 ||
+ param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
+ double rij,rik,costheta,lp1,lp3,lp6;
+ double rmu,rmu2,comtt,fck;
+ double pplp1 = param->plp1, pplp3 = param->plp3, pplp6 = param->plp6;
+ double c123 = cos(param->a123*PI/180.0);
+
+ // cos(theta) of the i-j-k
+ // cutoff function of rik
+
+ rij = sqrt(rsqij);
+ rik = sqrt(rsqik);
+ costheta = vec3_dot(delrij,delrik) / (rij*rik);
+ fck = comb_fc(rik,param);
+ rmu = costheta;
+
+ // Legendre Polynomial functions
+
+ if (param->plp1 > 1.0e-6 || param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
+ rmu2 = rmu*rmu;
+ lp1 = rmu; lp3 = 0.5*(5.0*rmu2*rmu-3.0*rmu);
+ lp6 = (231.0*rmu2*rmu2*rmu2-315.0*rmu2*rmu2+105.0*rmu2-5.0)/16.0;
+ comtt = pplp1*lp1 + pplp3*lp3 + pplp6*lp6;
+ } else comtt = 0.0;
+
+ // bond-bending terms
+
+ if (param->aconf>1e-4) {
+ if (param->hfocor >= 0.0)
+ comtt += param->aconf *(rmu-c123)*(rmu-c123);
+ else if (param->hfocor < 0.0)
+ comtt += param->aconf *(4.0-(rmu-c123)*(rmu-c123));
+ }
+
+ return 0.5 * fck * comtt;
+ }
+
+ return 0.0;
+}
+
+/* ----------------------------------------------------------------------
+ Legendre polynomial bond angle correction to forces
+------------------------------------------------------------------------- */
+
+void PairComb::flp(Param *param, double rsqij, double rsqik,
+ double *delrij, double *delrik, double *drilp,
+ double *drjlp, double *drklp)
+{
+ double ffj1,ffj2,ffk1,ffk2;
+ ffj1 = 0.0; ffj2 = 0.0; ffk1 = 0.0; ffk2 = 0.0;
+
+ if (param->aconf > 1.0e-4 || param->plp1 > 1.0e-6 ||
+ param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
+ double rij,rik,costheta,lp1,lp1_d,lp3,lp3_d,lp6,lp6_d;
+ double rmu,rmu2,comtt,comtt_d,com4k,com5,fcj,fck,fck_d;
+
+ double pplp1 = param->plp1;
+ double pplp3 = param->plp3;
+ double pplp6 = param->plp6;
+ double c123 = cos(param->a123*PI/180.0);
+
+ // fck_d = derivative of cutoff function
+
+ rij = sqrt(rsqij); rik = sqrt(rsqik);
+ costheta = vec3_dot(delrij,delrik) / (rij*rik);
+ fcj = comb_fc(rij,param);
+ fck = comb_fc(rik,param);
+ fck_d = comb_fc_d(rik,param);
+ rmu = costheta;
+
+ // Legendre Polynomial functions and derivatives
+
+ if (param->plp1 > 1.0e-6 || param->plp3 > 1.0e-6 || param->plp6 > 1.0e-6) {
+ rmu2 = rmu*rmu;
+ lp1 = rmu; lp3 = (2.5*rmu2*rmu-1.5*rmu);
+ lp6 = (231.0*rmu2*rmu2*rmu2-315.0*rmu2*rmu2+105.0*rmu2-5.0)/16.0;
+ lp1_d = 1.0;lp3_d = (7.5*rmu2-1.5);
+ lp6_d = (1386.0*rmu2*rmu2*rmu-1260.0*rmu2*rmu+210.0)/16.0;
+ comtt = pplp1*lp1 + pplp3*lp3 + pplp6*lp6;
+ comtt_d = pplp1*lp1_d + pplp3*lp3_d + pplp6*lp6_d;
+ } else {
+ comtt = 0.0;
+ comtt_d = 0.0;
+ }
+
+ // bond-bending terms derivatives
+
+ if (param->aconf > 1.0e-4) {
+ if (param->hfocor >= 0.0) {
+ comtt += param->aconf *(rmu-c123)*(rmu-c123);
+ comtt_d += 2.0*param->aconf*(rmu-c123);
+ } else if (param->hfocor < 0.0) {
+ comtt += param->aconf *(4.0-(rmu-c123)*(rmu-c123));
+ comtt_d += -2.0*param->aconf*(rmu-c123);
+ }
+ }
+
+ com4k = fcj * fck_d * comtt;
+ com5 = fcj * fck * comtt_d;
+
+ ffj1 =-0.5*(com5/(rij*rik));
+ ffj2 = 0.5*(com5*rmu/rsqij);
+ ffk1 = ffj1;
+ ffk2 = 0.5*(-com4k/rik+com5*rmu/rsqik);
+
+ } else {
+ ffj1 = 0.0; ffj2 = 0.0;
+ ffk1 = 0.0; ffk2 = 0.0;
+ }
+
+ // j-atom
+
+ vec3_scale(ffj1,delrik,drjlp); // (k,x[],y[]), y[]=k*x[]
+ vec3_scaleadd(ffj2,delrij,drjlp,drjlp); // (k,x[],y[],z[]), z[]=k*x[]+y[]
+
+ // k-atom
+
+ vec3_scale(ffk1,delrij,drklp);
+ vec3_scaleadd(ffk2,delrik,drklp,drklp);
+
+ // i-atom
+
+ vec3_add(drjlp,drklp,drilp); // (x[],y[],z[]), z[]=x[]+y[]
+ vec3_scale(-1.0,drilp,drilp);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::force_zeta(Param *param, double rsq, double zeta_ij,
+ double &fforce, double &prefactor,
+ int eflag, double &eng, double iq, double jq)
+{
+ double r,fa,fa_d,bij;
+
+ r = sqrt(rsq);
+ if (r > param->bigr+param->bigd) return;
+ fa = comb_fa(r,param,iq,jq);
+ fa_d = comb_fa_d(r,param,iq,jq);
+ bij = comb_bij(zeta_ij,param);
+ fforce = 0.5*bij*fa_d / r;
+ prefactor = -0.5*fa * comb_bij_d(zeta_ij,param);
+
+ // eng = attractive energy
+
+ if (eflag) eng = 0.5*bij*fa;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fc(double r, Param *param)
+{
+ double comb_R = param->bigr;
+ double comb_D = param->bigd;
+
+ if (r < comb_R-comb_D) return 1.0;
+ if (r > comb_R+comb_D) return 0.0;
+ return 0.5*(1.0 + cos(PI*(r - comb_R)/comb_D));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fc_d(double r, Param *param)
+{
+ double comb_R = param->bigr;
+ double comb_D = param->bigd;
+
+ if (r < comb_R-comb_D) return 0.0;
+ if (r > comb_R+comb_D) return 0.0;
+ return -(PI2/comb_D) * sin(PI*(r - comb_R)/comb_D);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fc2(double r)
+{
+ double comb_R = 1.89350;
+ double comb_D = comb_R + 0.050;
+
+ if (r < comb_R) return 0.0;
+ if (r > comb_D) return 1.0;
+ return 0.5*(1.0 + cos(PI*(r - comb_R)/(comb_D-comb_R)));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fc2_d(double r)
+{
+ double comb_R = 1.89350;
+ double comb_D = comb_R + 0.050;
+
+ if (r < comb_R) return 0.0;
+ if (r > comb_D) return 0.0;
+ return -(PI2/(comb_D-comb_R)) * sin(PI*(r - comb_R)/(comb_D-comb_R));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fc3(double r)
+{
+ double comb_R = 2.51350;
+ double comb_D = comb_R + 0.050;
+
+ if (r < comb_R) return 1.0;
+ if (r > comb_D) return 0.0;
+ return 0.5*(1.0 + cos(PI*(r - comb_R)/(comb_D-comb_R)));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fc3_d(double r)
+{
+ double comb_R = 2.51350;
+ double comb_D = comb_R + 0.050;
+
+ if (r < comb_R) return 0.0;
+ if (r > comb_D) return 0.0;
+ return -(PI2/(comb_D-comb_R)) * sin(PI*(r - comb_R)/(comb_D-comb_R));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::self(Param *param, double qi, double selfpot)
+{
+ double self_tmp, cmin, cmax, qmin, qmax;
+ double s1=param->chi, s2=param->dj, s3=param->dk, s4=param->dl, s5=param->dm;
+
+ self_tmp = 0.0;
+ qmin = param->QL1*0.90;
+ qmax = param->QU1*0.90;
+ cmin = cmax = 1000.0;
+
+ self_tmp = qi*(s1+qi*(s2+selfpot+qi*(s3+qi*(s4+qi*qi*s5))));
+
+ if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4);
+ if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4);
+
+ return self_tmp;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fa(double r, Param *param, double iq, double jq)
+{
+ double bigB,Bsi,Bsj;
+ double qi,qj,Di,Dj;
+
+ if (r > param->bigr + param->bigd) return 0.0;
+ qi = iq; qj = jq;
+ Di = Dj = Bsi = Bsj = bigB = 0.0;
+ Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1);
+ Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2);
+ Bsi = param->bigb1 * exp(param->lam21*Di)*
+ (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10)));
+ Bsj = param->bigb2 * exp(param->lam22*Dj)*
+ (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10)));
+ if (Bsi > 0.0 && Bsj > 0.0) bigB = sqrt(Bsi*Bsj)*param->romigb;
+ else bigB = 0.0;
+
+ return -bigB * exp(-param->rlm2 * r) * comb_fc(r,param);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_fa_d(double r, Param *param, double iq, double jq)
+{
+ double bigB,Bsi,Bsj;
+ double qi,qj,Di,Dj;
+
+ if (r > param->bigr + param->bigd) return 0.0;
+ qi = iq; qj = jq;
+ Di = Dj = Bsi = Bsj = bigB = 0.0;
+ Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1);
+ Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2);
+ Bsi = param->bigb1 * exp(param->lam21*Di)*
+ (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10)));
+ Bsj = param->bigb2 * exp(param->lam22*Dj)*
+ (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10)));
+ if (Bsi > 0.0 && Bsj > 0.0) bigB = sqrt(Bsi*Bsj)*param->romigb;
+ else bigB = 0.0;
+
+ return bigB * exp(-param->rlm2 * r) *
+ (param->rlm2 * comb_fc(r,param) - comb_fc_d(r,param));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_bij(double zeta, Param *param)
+{
+ double tmp = param->beta * zeta;
+ if (tmp > param->c1) return 1.0/sqrt(tmp);
+ if (tmp > param->c2)
+ return (1.0 - pow(tmp,-param->powern) / (2.0*param->powern))/sqrt(tmp);
+ if (tmp < param->c4) return 1.0;
+ if (tmp < param->c3)
+ return 1.0 - pow(tmp,param->powern)/(2.0*param->powern);
+ return pow(1.0 + pow(tmp,param->powern), -1.0/(2.0*param->powern));
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_bij_d(double zeta, Param *param)
+{
+ double tmp = param->beta * zeta;
+ if (tmp > param->c1) return param->beta * -0.5*pow(tmp,-1.5);
+ if (tmp > param->c2)
+ return param->beta * (-0.5*pow(tmp,-1.5) *
+ (1.0 - 0.5*(1.0 + 1.0/(2.0*param->powern)) *
+ pow(tmp,-param->powern)));
+ if (tmp < param->c4) return 0.0;
+ if (tmp < param->c3)
+ return -0.5*param->beta * pow(tmp,param->powern-1.0);
+
+ double tmp_n = pow(tmp,param->powern);
+ return -0.5 * pow(1.0+tmp_n, -1.0-(1.0/(2.0*param->powern)))*tmp_n / zeta;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_gijk(double costheta, Param *param)
+{
+ double comb_c = param->c;
+ double comb_d = param->d;
+
+ return (1.0 + pow(comb_c/comb_d,2.0) -
+ pow(comb_c,2.0) / (pow(comb_d,2.0) + pow(param->h - costheta,2.0)));
+};
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::comb_gijk_d(double costheta, Param *param)
+{
+ double numerator = -2.0 * pow(param->c,2) * (param->h - costheta);
+ double denominator = pow(pow(param->d,2.0) +
+ pow(param->h - costheta,2.0),2.0);
+ return numerator/denominator;
+}
+
+/*------------------------------------------------------------------------- */
+
+void PairComb::attractive(Param *param, double prefactor,
+ double rsqij, double rsqik,
+ double *delrij, double *delrik,
+ double *fi, double *fj, double *fk)
+{
+ double rij_hat[3],rik_hat[3];
+ double rij,rijinv,rik,rikinv;
+
+ rij = sqrt(rsqij);
+ rijinv = 1.0/rij;
+ vec3_scale(rijinv,delrij,rij_hat);
+
+ rik = sqrt(rsqik);
+ rikinv = 1.0/rik;
+ vec3_scale(rikinv,delrik,rik_hat);
+
+ comb_zetaterm_d(prefactor,rij_hat,rij,rik_hat,rik,fi,fj,fk,param);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::comb_zetaterm_d(double prefactor, double *rij_hat, double rij,
+ double *rik_hat, double rik, double *dri,
+ double *drj, double *drk, Param *param)
+{
+ double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta,tmp;
+ double dcosdri[3],dcosdrj[3],dcosdrk[3];
+
+ fc = comb_fc(rik,param);
+ dfc = comb_fc_d(rik,param);
+ if (param->powermint == 3) tmp = pow(param->rlm2 * (rij-rik),3.0);
+ else tmp = param->rlm2 * (rij-rik);
+
+ if (tmp > 69.0776) ex_delr = 1.e30;
+ else if (tmp < -69.0776) ex_delr = 0.0;
+ else ex_delr = exp(tmp); // ex_delr is Ygexp
+
+ if (param->powermint == 3)
+ ex_delr_d = 3.0*pow(param->rlm2,3.0) * pow(rij-rik,2.0)*ex_delr; // com3
+ else ex_delr_d = param->rlm2 * ex_delr; // com3
+
+ cos_theta = vec3_dot(rij_hat,rik_hat);
+ gijk = comb_gijk(cos_theta,param);
+ gijk_d = comb_gijk_d(cos_theta,param);
+ 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);
+ // (k,x[],y[]), y[]=k*x[]
+ // (k,x[],y[],z[]), z[]=k*x[]+y[]
+
+ 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 PairComb::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);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::sm_table()
+{
+ int i,j,k,m,nntypes,ncoul;
+ int inty, itype, jtype;
+ int iparam_i, iparam_ij, iparam_ji;
+ double r,dra,drin,rc,z,zr,zrc,ea,eb,ea3,eb3,alf;
+ double exp2er,exp2ersh,fafash,dfafash,F1,dF1,ddF1,E1,E2,E3,E4;
+ double exp2ear,exp2ebr,exp2earsh,exp2ebrsh,fafbsh,dfafbsh;
+
+ int n = atom->ntypes;
+ int *type = atom->type;
+
+ dra = 0.001; // lookup table step size
+ drin = 0.1; // starting distance of 1/r
+ rc = cutmax;
+ alf = 0.20;
+
+ nntypes = int((n+1)*n/2); // interaction types
+ ncoul = int((rc-drin)/dra)+1;
+
+ memory->destroy_2d_int_array(intype);
+ memory->destroy_2d_double_array(fafb);
+ memory->destroy_2d_double_array(dfafb);
+ memory->destroy_2d_double_array(ddfafb);
+ memory->destroy_2d_double_array(phin);
+ memory->destroy_2d_double_array(dphin);
+ memory->destroy_2d_double_array(erpaw);
+
+ // allocate arrays
+
+ intype = memory->create_2d_int_array(n,n,"pair:intype");
+ fafb = memory->create_2d_double_array(ncoul,nntypes,"pair:fafb");
+ dfafb = memory->create_2d_double_array(ncoul,nntypes,"pair:dfafb");
+ ddfafb = memory->create_2d_double_array(ncoul,nntypes,"pair:ddfafb");
+ phin = memory->create_2d_double_array(ncoul,nntypes,"pair:phin");
+ dphin = memory->create_2d_double_array(ncoul,nntypes,"pair:dphin");
+ erpaw = memory->create_2d_double_array(25000,2,"pair:erpaw");
+
+ // set interaction number: 0-0=0, 1-1=1, 0-1=1-0=2
+
+ m = 0; k = n;
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < n; j++) {
+ if (j == i) {
+ intype[i][j] = m;
+ m += 1;
+ } else if (j != i && j > i) {
+ intype[i][j] = k;
+ k += 1;
+ } else if (j != i && j < i) {
+ intype[i][j] = intype[j][i];
+ }
+ }
+ }
+
+ // default arrays to zero
+
+ for (i = 0; i < ncoul; i ++) {
+ for (j = 0; j < nntypes; j ++) {
+ fafb[i][j] = 0.0;
+ dfafb[i][j] = 0.0;
+ ddfafb[i][j] = 0.0;
+ phin[i][j] = 0.0;
+ dphin[i][j] = 0.0;
+ }
+ }
+
+ // direct 1/r energy with Slater 1S orbital overlap
+
+ for (i = 0; i < n; i++) {
+ r = drin;
+ itype = params[i].ielement;
+ iparam_i = elem2param[itype][itype][itype];
+ z = params[iparam_i].esm1;
+ for (j = 0; j < ncoul; j++) {
+ exp2er = exp(-2.0 * z * r);
+ phin[j][i] = 1.0 - exp2er * (1.0 + 2.0 * z * r * (1.0 + z * r));
+ dphin[j][i] = (4.0 * exp2er * z * z * z * r * r);
+ r += dra;
+ }
+ }
+
+ for (i = 0; i < n; i ++) {
+ for (j = 0; j < n; j ++) {
+ r = drin;
+ if (j == i) {
+ itype = params[i].ielement;
+ inty = intype[itype][itype];
+ iparam_i = elem2param[itype][itype][itype];
+ z = params[iparam_i].esm1;
+ zrc = z * rc;
+ exp2ersh = exp(-2.0 * zrc);
+ fafash = -exp2ersh * (1.0 / rc +
+ z * (11.0/8.0 + 3.0/4.0*zrc + zrc*zrc/6.0));
+ dfafash = exp2ersh * (1.0/(rc*rc) + 2.0*z/rc +
+ z*z*(2.0 + 7.0/6.0*zrc + zrc*zrc/3.0));
+ for (k = 0; k < ncoul; k ++) {
+ zr = z * r;
+ exp2er = exp(-2.0*zr);
+ F1 = -exp2er * (1.0 / r +
+ z * (11.0/8.0 + 3.0/4.0*zr + zr*zr/6.0));
+ dF1 = exp2er * (1.0/(r*r) + 2.0*z/r +
+ z*z*(2.0 + 7.0/6.0*zr + zr*zr/3.0));
+ ddF1 = -exp2er * (2.0/(r*r*r) + 4.0*z/(r*r) -
+ z*z*z/3.0*(17.0/2.0 + 5.0*zr + 2.0*zr*zr));
+ fafb[k][inty] = F1-fafash-(r-rc)*dfafash;
+ dfafb[k][inty] = (dF1 - dfafash);
+ ddfafb[k][inty] = ddF1;
+ r += dra;
+ }
+ } else if (j != i) {
+ itype = params[i].ielement;
+ jtype = params[j].ielement;
+ inty = intype[itype][jtype];
+ iparam_ij = elem2param[itype][jtype][jtype];
+ ea = params[iparam_ij].esm1;
+ ea3 = ea*ea*ea;
+ iparam_ji = elem2param[jtype][itype][itype];
+ eb = params[iparam_ji].esm1;
+ eb3 = eb*eb*eb;
+ E1 = ea*eb3*eb/((ea+eb)*(ea+eb)*(ea-eb)*(ea-eb));
+ E2 = eb*ea3*ea/((ea+eb)*(ea+eb)*(eb-ea)*(eb-ea));
+ E3 = (3.0*ea*ea*eb3*eb-eb3*eb3) /
+ ((ea+eb)*(ea+eb)*(ea+eb)*(ea-eb)*(ea-eb)*(ea-eb));
+ E4 = (3.0*eb*eb*ea3*ea-ea3*ea3) /
+ ((ea+eb)*(ea+eb)*(ea+eb)*(eb-ea)*(eb-ea)*(eb-ea));
+ exp2earsh = exp(-2.0*ea*rc);
+ exp2ebrsh = exp(-2.0*eb*rc);
+ fafbsh = -exp2earsh*(E1 + E3/rc)-exp2ebrsh*(E2 + E4/rc);
+ dfafbsh =
+ exp2earsh*(2.0*ea*(E1+E3/rc)+E3/(rc*rc)) +
+ exp2ebrsh*(2.0*eb*(E2+E4/rc)+E4/(rc*rc));
+ for (k = 0; k < ncoul; k ++) {
+ exp2ear = exp(-2.0*ea*r);
+ exp2ebr = exp(-2.0*eb*r);
+ fafb[k][inty] =
+ - exp2ear*(E1+E3/r) - exp2ebr*(E2+E4/r)
+ - fafbsh - (r-rc) * dfafbsh;
+ dfafb[k][inty] = (exp2ear*(2.0*ea*(E1+E3/r) + E3/(r*r))
+ + exp2ebr*(2.0*eb*(E2+E4/r) + E4/(r*r))- dfafbsh);
+ ddfafb[k][inty] = (- exp2ear*(E3/(r*r)*(1.0/r+2.0*ea/r+2.0/(r*r))
+ + 2.0*ea*(E1+E3/r))-
+ exp2ebr*(E4/(r*r)
+ *(1.0/r+2.0*eb/r+2.0/(r*r)) +
+ 2.0*eb*(E2+E4/r)));
+ r += dra;
+ }
+ }
+ }
+ }
+
+ for (i = 0; i < 25000; i ++) {
+ r = dra * i + drin;
+ erpaw[i][0] = erfc(r*alf);
+ erpaw[i][1] = exp(-r*r*alf*alf);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::potal_calc(double &calc1, double &calc2, double &calc3)
+{
+ double potal,alf,rcoul,fac11,fac11e,esucon;
+ int m;
+
+ rcoul = 0.0;
+ for (m = 0; m < nparams; m++)
+ if (params[m].lcut > rcoul) rcoul = params[m].lcut;
+
+ alf = 0.20;
+ esucon = force->qqr2e;
+
+ calc2 = (erfc(rcoul*alf)/rcoul/rcoul+2.0*alf/PIsq*
+ exp(-alf*alf*rcoul*rcoul)/rcoul)*esucon/rcoul;
+ calc3 = (erfc(rcoul*alf)/rcoul)*esucon;
+ calc1 = -(alf/PIsq*esucon+calc3*0.5);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::tri_point(double rsq, int &mr1, int &mr2,
+ int &mr3, double &sr1, double &sr2,
+ double &sr3, int &itype)
+{
+ double r, rin, dr, dd, rr1, rridr, rridr2;
+ int m = itype;
+
+ rin = 0.10; dr = 0.0010;
+ r = sqrt(rsq);
+ if (r < rin + 2.0*dr) r = rin + 2.0*dr;
+ if (r > cutmax - 2.0*dr) r = cutmax - 2.0*dr;
+ rridr = (r-rin)/dr;
+
+ mr1 = int(rridr)-1;
+ dd = rridr - float(mr1);
+ if (dd > 0.5) mr1 += 1;
+ mr2 = mr1 + 1;
+ mr3 = mr2 + 1;
+
+ rr1 = float(mr1)*dr;
+ rridr = (r - rin - rr1)/dr;
+ rridr2 = rridr * rridr;
+
+ sr1 = (rridr2 - rridr) * 0.50;
+ sr2 = 1.0 - rridr2;
+ sr3 = (rridr2 + rridr) * 0.50;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::direct(int inty, int mr1, int mr2, int mr3, double rsq,
+ double sr1, double sr2, double sr3,
+ double iq, double jq,
+ double potal, double fac11, double fac11e,
+ double &pot_tmp, double &pot_d)
+{
+ double r,erfcc,fafbn1,potij,sme2,chrij,esucon;
+ double r3,erfcd,dfafbn1,smf2,dvdrr,alf,alfdpi;
+
+ r = sqrt(rsq);
+ r3 = r * rsq;
+ alf = 0.20;
+ alfdpi = 2.0*alf/PIsq;
+ esucon = force->qqr2e;
+ pot_tmp = 0.0;
+ pot_d = 0.0;
+
+ // 1/r energy
+
+ erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
+ fafbn1= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
+ potij = (erfcc/r * esucon - fac11e);
+ sme2 = potij + fafbn1 * esucon;
+ pot_tmp = 1.0 * iq * jq *sme2;
+
+ // 1/r force (wrt r)
+
+ erfcd = sr1*erpaw[mr1][1] + sr2*erpaw[mr2][1] + sr3*erpaw[mr3][1];
+ dfafbn1= sr1*dfafb[mr1][inty] + sr2*dfafb[mr2][inty] + sr3*dfafb[mr3][inty];
+ dvdrr = (erfcc/r3+alfdpi*erfcd/rsq)*esucon-fac11;
+ smf2 = dvdrr - dfafbn1 * esucon/r;
+ pot_d = 1.0 * iq * jq * smf2;
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::field(Param *param, double rsq, double iq,double jq,
+ double &vionij,double &fvionij)
+{
+ double r,r5,r6,rc,rc5,rc6,rf5,drf6,smpn,smpl,rfx1,rfx2;
+ double cmi1,cmi2,cmj1,cmj2;
+
+ r = sqrt(rsq);
+ r5 = r*r*r*r*r;
+ r6 = r5 * r;
+ rc = param->lcut;
+ rc5 = rc*rc*rc*rc*rc;
+ rc6 = rc5 * rc;
+ cmi1 = param->cmn1;
+ cmi2 = param->cmn2;
+ cmj1 = param->cml1;
+ cmj2 = param->cml2;
+ rf5 = 1.0/r5 - 1.0/rc5 + 5.0*(r-rc)/rc6;
+ drf6 = 5.0/rc6 - 5.0/r6;
+
+ // field correction energy
+
+ smpn = rf5*jq*(cmi1+jq*cmi2);
+ smpl = rf5*iq*(cmj1+iq*cmj2);
+ vionij += 1.0 * (smpn + smpl);
+
+ // field correction force
+
+ rfx1 = jq*drf6*(cmi1+jq*cmi2)/r;
+ rfx2 = iq*drf6*(cmj1+iq*cmj2)/r;
+ fvionij -= 1.0 * (rfx1 + rfx2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::yasu_char(double *qf_fix, int &igroup)
+{
+ int i,j,k,ii,jj,kk,jnum;
+ int itag,jtag,itype,jtype,ktype,iparam_i,iparam_ij,iparam_ijk;
+ double xtmp,ytmp,ztmp,evdwl,fpair;
+ double rsq,rsq1,rsq2,delr1[3],delr2[3],zeta_ij;
+ int *ilist,*jlist,*numneigh,**firstneigh;
+ double iq,jq,fqi,fqj,fqij,fqjj,ecoul,yaself,yaself_d;
+ double potal,fac11,fac11e,sr1,sr2,sr3;
+ int mr1,mr2,mr3,inty;
+ int zeta_flag;
+
+ double **x = atom->x;
+ double *q = atom->q;
+ int *tag = atom->tag;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ int inum = list->inum;
+ ilist = list->ilist;
+ numneigh = list->numneigh;
+ firstneigh = list->firstneigh;
+
+ int *mask = atom->mask;
+ int groupbit = group->bitmask[igroup];
+
+ qf = qf_fix;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (mask[i] & groupbit)
+ qf[i] = 0.0;
+ }
+
+ // communicating charge force to all nodes, first forward then reverse
+
+ comm->forward_comm_pair(this);
+
+ // self energy correction term: potal
+
+ potal_calc(potal,fac11,fac11e);
+
+ // loop over full neighbor list of my atoms
+
+ fqi = fqj = fqij = fqjj = 0.0;
+
+ for (ii = 0; ii < inum; ii ++) {
+ i = ilist[ii];
+ if (mask[i] & groupbit) {
+ itype = map[type[i]];
+ xtmp = x[i][0];
+ ytmp = x[i][1];
+ ztmp = x[i][2];
+ iq = q[i];
+ iparam_i = elem2param[itype][itype][itype];
+
+ // charge force from self energy
+
+ fqi = qfo_self(¶ms[iparam_i],iq,potal);
+
+ // two-body interactions
+
+ jlist = firstneigh[i];
+ jnum = numneigh[i];
+
+ for (jj = 0; jj < jnum; jj++) {
+ j = jlist[jj];
+ jtype = map[type[j]];
+ jq = q[j];
+
+ delr1[0] = x[j][0] - xtmp;
+ delr1[1] = x[j][1] - ytmp;
+ delr1[2] = x[j][2] - ztmp;
+ rsq1 = vec3_dot(delr1,delr1);
+
+ iparam_ij = elem2param[itype][jtype][jtype];
+
+ // long range q-dependent
+
+ if (rsq1 > params[iparam_ij].lcutsq) continue;
+
+ inty = intype[itype][jtype];
+
+ // polynomial three-point interpolation
+
+ tri_point(rsq1,mr1,mr2,mr3,sr1,sr2,sr3,itype);
+
+ // 1/r charge forces
+
+ qfo_direct(inty,mr1,mr2,mr3,rsq1,sr1,sr2,sr3,fac11e,fqij);
+ fqi += jq * fqij; qf[j] += iq * fqij;
+
+ // field correction to self energy and charge force
+
+ qfo_field(¶ms[iparam_ij],rsq1,iq,jq,fqij,fqjj);
+ fqi += fqij;
+ qf[j] += fqjj;
+
+ // polarization field charge force
+ // three-body interactions
+
+ if (rsq1 > params[iparam_ij].cutsq) continue;
+
+ zeta_ij = 0.0;
+
+ for (kk = 0; kk < jnum; kk++) {
+ if (jj == kk) continue;
+ k = jlist[kk];
+ ktype = map[type[k]];
+ iparam_ijk = elem2param[itype][jtype][ktype];
+
+ delr2[0] = x[k][0] - xtmp;
+ delr2[1] = x[k][1] - ytmp;
+ delr2[2] = x[k][2] - ztmp;
+ rsq2 = vec3_dot(delr2,delr2);
+
+ if (rsq2 > params[iparam_ijk].cutsq) continue;
+ zeta_ij += zeta(¶ms[iparam_ijk],rsq1,rsq2,delr1,delr2);
+ }
+
+ // charge force in Aij and Bij
+
+ qfo_short(¶ms[iparam_ij],rsq1,zeta_ij,iq,jq,fqij,fqjj);
+ fqi += fqij; qf[j] += fqjj;
+ }
+
+ qf[i] += fqi;
+
+ }
+ }
+
+ comm->reverse_comm_pair(this);
+
+ // sum charge force on each node and return it
+
+ double eneg = 0.0;
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (mask[i] & groupbit)
+ eneg += qf[i];
+ }
+ double enegtot;
+ MPI_Allreduce(&eneg,&enegtot,1,MPI_DOUBLE,MPI_SUM,world);
+ return enegtot;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double PairComb::qfo_self(Param *param, double qi, double selfpot)
+{
+ double self_d,cmin,cmax,qmin,qmax;
+ double s1 = param->chi;
+ double s2 = param->dj;
+ double s3 = param->dk;
+ double s4 = param->dl;
+ double s5 = param->dm;
+
+ self_d = 0.0;
+ qmin = param->QL1*0.90;
+ qmax = param->QU1*0.90;
+ if (qmax > 4.50 ) qmax = -0.70;
+ cmin = cmax = 1000.0;
+
+ self_d = s1+qi*(2.0*(s2+selfpot)+qi*(3.0*s3+qi*(4.0*s4+qi*qi*6.0*s5)));
+
+ if (qi < qmin) {
+ // char str[128];
+ // sprintf(str,"Pair COMB charge %.10f with force %.10f hit min barrier",
+ // qi,self_d);
+ // error->warning(str,0);
+ self_d += 4.0 * cmin * pow((qi-qmin),3);
+ }
+ if (qi > qmax) {
+ // char str[128];
+ // sprintf(str,"Pair COMB charge %.10f with force %.10f hit max barrier",
+ // qi,self_d);
+ // error->warning(str,0);
+ self_d += 4.0 * cmax * pow((qi-qmax),3);
+ }
+
+ return self_d;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::qfo_direct(int inty, int mr1, int mr2, int mr3,
+ double rsq, double sr1, double sr2,
+ double sr3, double fac11e, double &fqij)
+{
+ double r, erfcc, fafbn1, vm, esucon;
+
+ r = sqrt(rsq);
+ esucon=force->qqr2e;
+
+ // 1/r force (wrt q)
+
+ erfcc = sr1*erpaw[mr1][0] + sr2*erpaw[mr2][0] + sr3*erpaw[mr3][0];
+ fafbn1= sr1*fafb[mr1][inty] + sr2*fafb[mr2][inty] + sr3*fafb[mr3][inty];
+ vm = (erfcc/r * esucon - fac11e);
+ fqij = 0.5 * (vm+esucon*fafbn1);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::qfo_field(Param *param, double rsq,double iq,double jq,
+ double &fqij, double &fqjj)
+{
+ double r,r5,r6,rc,rc5,rc6;
+ double cmi1,cmi2,cmj1,cmj2,rf5;
+
+ fqij = fqjj = 0.0;
+ r = sqrt(rsq);
+ r5 = r*r*r*r*r;
+ r6 = r5 * r;
+ rc = param->lcut;
+ rc5 = rc*rc*rc*rc*rc;
+ rc6 = rc5 * rc;
+ cmi1 = param->cmn1;
+ cmi2 = param->cmn2;
+ cmj1 = param->cml1;
+ cmj2 = param->cml2;
+ rf5 = 1.0/r5 - 1.0/rc5 + 5.0*(r-rc)/rc6;
+
+ // field correction charge force
+
+ fqij = 0.5 * rf5 * (cmj1 + 2.0 * iq * cmj2);
+ fqjj = 0.5 * rf5 * (cmi1 + 2.0 * jq * cmi2);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::qfo_short(Param *param, double rsq, double zeta_ij,
+ double iq, double jq, double &fqij, double &fqjj)
+{
+ double r,tmp_fc,tmp_fc_d,tmp_exp1,tmp_exp2;
+ double bigA,Asi,Asj,vrcs;
+ double romi = param->addrep,rrcs = param->bigr + param->bigd;
+ double qi,qj,Di,Dj,bigB,Bsi,Bsj;
+ double QUchi,QOchi,QUchj,QOchj,YYDiqp,YYDjqp;
+ double YYAsiqp,YYAsjqp,YYBsiqp,YYBsjqp;
+ double caj,cbj,bij,cfqr,cfqs;
+ double romie = param->romiga;
+ double romib = param->romigb;
+ double ca1,ca2,ca3,ca4;
+ double rslp,rslp2,rslp4,arr1,arr2,fc2j,fc3j,fcp2j,fcp3j;
+
+ qi = iq; qj = jq; r = sqrt(rsq);
+ Di = Dj = Asi = Asj = bigA = Bsi = Bsj = bigB = 0.0;
+ QUchi = QOchi = QUchj = QOchj = YYDiqp = YYDjqp =0.0;
+ YYAsiqp = YYAsjqp = YYBsiqp = YYBsjqp = 0.0;
+ caj = cbj = vrcs = cfqr = cfqs = 0.0;
+
+ tmp_fc = comb_fc(r,param);
+ tmp_fc_d = comb_fc_d(r,param);
+ tmp_exp1 = exp(-param->rlm1 * r);
+ tmp_exp2 = exp(-param->rlm2 * r);
+ bij = comb_bij(zeta_ij,param);
+
+ arr1 = 2.22850; arr2 = 1.89350;
+ fc2j = comb_fc2(r);
+ fc3j = comb_fc3(r);
+ fcp2j = comb_fc2_d(r);
+ fcp3j = comb_fc3_d(r);
+
+ vrcs = 0.0;
+ if (romi > 0.0) {
+ if (!cor_flag) vrcs = romi * pow((1.0-r/rrcs),2.0);
+ else if (cor_flag) {
+ rslp = ((arr1-r)/(arr1-arr2));
+ rslp2 = rslp * rslp; rslp4 = rslp2 * rslp2;
+ vrcs = fc2j * fc3j * romi * ((50.0*rslp4-30.0*rslp2+4.50))/8.0;
+ }
+ }
+
+ Di = param->DU1 + pow(fabs(param->bD1*(param->QU1-qi)),param->nD1);
+ Dj = param->DU2 + pow(fabs(param->bD2*(param->QU2-qj)),param->nD2);
+
+ Asi = param->biga1 * exp(param->lam11*Di);
+ Asj = param->biga2 * exp(param->lam12*Dj);
+ Bsi = param->bigb1 * exp(param->lam21*Di)*
+ (param->aB1-fabs(pow(param->bB1*(qi-param->Qo1),10)));
+ Bsj = param->bigb2 * exp(param->lam22*Dj)*
+ (param->aB2-fabs(pow(param->bB2*(qj-param->Qo2),10)));
+
+ QUchi = (param->QU1-qi)*param->bD1;
+ QUchj = (param->QU2-qj)*param->bD2;
+ QOchi = (qi-param->Qo1)*param->bB1;
+ QOchj = (qj-param->Qo2)*param->bB2;
+
+ if (QUchi == 0.0) YYDiqp = 0.0;
+ else YYDiqp = -param->nD1 * QUchi * param->bD1 *
+ pow(fabs(QUchi),(param->nD1-2.0));
+
+ if (QUchj == 0.0) YYDjqp = 0.0;
+ else YYDjqp = -param->nD2 * QUchj * param->bD2 *
+ pow(fabs(QUchj),(param->nD2-2.0));
+
+ YYAsiqp = Asi * param->lam11 * YYDiqp;
+ YYAsjqp = Asj * param->lam12 * YYDjqp;
+
+ if (QOchi == 0.0)
+ YYBsiqp=Bsi*param->lam21*YYDiqp;
+ else
+ YYBsiqp=Bsi*param->lam21*YYDiqp-param->bigb1*exp(param->lam21*Di)*
+ 10.0*QOchi*param->bB1*pow(fabs(QOchi),(10.0-2.0));
+
+ if (QOchj == 0.0)
+ YYBsjqp=Bsj*param->lam22*YYDjqp;
+ else
+ YYBsjqp=Bsj*param->lam22*YYDjqp-param->bigb2*exp(param->lam22*Dj)*
+ 10.0*QOchj*param->bB2*pow(fabs(QOchj),(10.0-2.0));
+
+ if (Asi > 0.0 && Asj > 0.0) caj = 1.0/(2.0*sqrt(Asi*Asj)) * romie;
+ else caj == 0.0;
+
+ if (Bsi > 0.0 && Bsj > 0.0) cbj = 1.0/(2.0*sqrt(Bsi*Bsj)) * romib ;
+ else cbj == 0.0;
+
+ cfqr = 0.50 * tmp_fc * (1.0 + vrcs); // 0.5 b/c full atom loop
+ cfqs = -0.50 * tmp_fc * bij;
+
+ ca1 = Asj * caj * YYAsiqp;
+ ca2 = Bsj * cbj * YYBsiqp;
+ ca3 = Asi * caj * YYAsjqp;
+ ca4 = Bsi * cbj * YYBsjqp;
+
+ fqij = cfqr * tmp_exp1 * ca1;
+ fqij += cfqs * tmp_exp2 * ca2;
+ fqjj = cfqr * tmp_exp1 * ca3;
+ fqjj += cfqs * tmp_exp2 * ca4;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::Over_cor(Param *param, double rsq1, int NCoi,
+ double &Eov, double &Fov)
+{
+ double ECo,BCo,tmp_fc,tmp_fc_d;
+ double r = sqrt(rsq1);
+ int NCon = NCoi - 7;
+
+ tmp_fc = comb_fc(r,param);
+ tmp_fc_d = comb_fc(r,param);
+ Eov = 0.0; Fov = 0.0;
+ ECo = param->hfocor;
+ BCo = 0.1;
+
+ if (NCon >= 0.20) {
+ Eov = tmp_fc * ECo * NCon/(1.0+exp(BCo*NCon));
+ Fov = -(tmp_fc_d*Eov + tmp_fc*ECo/(1.0+exp(BCo*NCon)) -
+ (tmp_fc*ECo*NCon*BCo*exp(BCo*NCon)) /
+ ((1.0+exp(BCo*NCon))*(1.0+exp(BCo*NCon))));
+ Fov /= r;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairComb::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i ++) {
+ j = list[i];
+ buf[m++] = qf[j];
+ }
+ return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::unpack_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n ;
+ for (i = first; i < last; i++) qf[i] = buf[m++];
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairComb::pack_reverse_comm(int n, int first, double *buf)
+{
+ int i,m,last;
+
+ m = 0;
+ last = first + n;
+ for (i = first; i < last; i++) buf[m++] = qf[i];
+ return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairComb::unpack_reverse_comm(int n, int *list, double *buf)
+{
+ int i,j,m;
+
+ m = 0;
+ for (i = 0; i < n; i++) {
+ j = list[i];
+ qf[j] += buf[m++];
+ }
+}
+
+/* ----------------------------------------------------------------------
+ memory usage of local atom-based arrays
+------------------------------------------------------------------------- */
+
+double PairComb::memory_usage()
+{
+ double bytes = maxeatom * sizeof(double);
+ bytes += maxvatom*6 * sizeof(double);
+ bytes += nmax * sizeof(int);
+ return bytes;
+}