diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp index 52d250fbe..63c0cb08f 100644 --- a/src/KSPACE/pair_buck_long_coul_long.cpp +++ b/src/KSPACE/pair_buck_long_coul_long.cpp @@ -1,1053 +1,1052 @@ /* ---------------------------------------------------------------------- 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: Pieter J. in 't Veld (SNL) ------------------------------------------------------------------------- */ #include #include #include #include #include "math_vector.h" #include "pair_buck_long_coul_long.h" #include "atom.h" #include "comm.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "force.h" #include "kspace.h" #include "update.h" #include "integrate.h" #include "respa.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define EWALD_F 1.12837917 #define EWALD_P 0.3275911 #define A1 0.254829592 #define A2 -0.284496736 #define A3 1.421413741 #define A4 -1.453152027 #define A5 1.061405429 /* ---------------------------------------------------------------------- */ PairBuckLongCoulLong::PairBuckLongCoulLong(LAMMPS *lmp) : Pair(lmp) { dispersionflag = ewaldflag = pppmflag = 1; respa_enable = 1; writedata = 1; - ftable = NULL; - fdisptable = NULL; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::options(char **arg, int order) { const char *option[] = {"long", "cut", "off", NULL}; int i; if (!*arg) error->all(FLERR,"Illegal pair_style buck/long/coul/long command"); for (i=0; option[i]&&strcmp(arg[0], option[i]); ++i); switch (i) { default: error->all(FLERR,"Illegal pair_style buck/long/coul/long command"); case 0: ewald_order |= 1<all(FLERR,"Illegal pair_style command"); ewald_order = 0; ewald_off = 0; options(arg,6); options(++arg,1); if (!comm->me && ewald_order == ((1<<1) | (1<<6))) error->warning(FLERR,"Using largest cutoff for buck/long/coul/long"); if (!*(++arg)) error->all(FLERR,"Cutoffs missing in pair_style buck/long/coul/long"); if (ewald_off & (1<<6)) error->all(FLERR,"LJ6 off not supported in pair_style buck/long/coul/long"); if (!((ewald_order^ewald_off) & (1<<1))) error->all(FLERR, "Coulomb cut not supported in pair_style buck/long/coul/coul"); cut_buck_global = force->numeric(FLERR,*(arg++)); if (narg == 4 && ((ewald_order & 0x42) == 0x42)) error->all(FLERR,"Only one cutoff allowed when requesting all long"); if (narg == 4) cut_coul = force->numeric(FLERR,*arg); else cut_coul = cut_buck_global; if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut_buck[i][j] = cut_buck_global; } } /* ---------------------------------------------------------------------- free all arrays ------------------------------------------------------------------------- */ PairBuckLongCoulLong::~PairBuckLongCoulLong() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut_buck_read); memory->destroy(cut_buck); memory->destroy(cut_bucksq); memory->destroy(buck_a_read); memory->destroy(buck_a); memory->destroy(buck_c_read); memory->destroy(buck_c); memory->destroy(buck_rho_read); memory->destroy(buck_rho); memory->destroy(buck1); memory->destroy(buck2); memory->destroy(rhoinv); memory->destroy(offset); } if (ftable) free_tables(); + if (fdisptable) free_disp_tables(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cut_buck_read,n+1,n+1,"pair:cut_buck_read"); memory->create(cut_buck,n+1,n+1,"pair:cut_buck"); memory->create(cut_bucksq,n+1,n+1,"pair:cut_bucksq"); memory->create(buck_a_read,n+1,n+1,"pair:buck_a_read"); memory->create(buck_a,n+1,n+1,"pair:buck_a"); memory->create(buck_c_read,n+1,n+1,"pair:buck_c_read"); memory->create(buck_c,n+1,n+1,"pair:buck_c"); memory->create(buck_rho_read,n+1,n+1,"pair:buck_rho_read"); memory->create(buck_rho,n+1,n+1,"pair:buck_rho"); memory->create(buck1,n+1,n+1,"pair:buck1"); memory->create(buck2,n+1,n+1,"pair:buck2"); memory->create(rhoinv,n+1,n+1,"pair:rhoinv"); memory->create(offset,n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- extract protected data from object ------------------------------------------------------------------------- */ void *PairBuckLongCoulLong::extract(const char *id, int &dim) { const char *ids[] = { "B", "ewald_order", "ewald_cut", "ewald_mix", "cut_coul", "cut_LJ", NULL}; void *ptrs[] = { buck_c, &ewald_order, &cut_coul, &mix_flag, &cut_coul, &cut_buck_global, NULL}; int i; for (i=0; ids[i]&&strcmp(ids[i], id); ++i); if (i == 0) dim = 2; else dim = 0; return ptrs[i]; } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::coeff(int narg, char **arg) { if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(*(arg++),atom->ntypes,ilo,ihi); force->bounds(*(arg++),atom->ntypes,jlo,jhi); double buck_a_one = force->numeric(FLERR,*(arg++)); double buck_rho_one = force->numeric(FLERR,*(arg++)); double buck_c_one = force->numeric(FLERR,*(arg++)); double cut_buck_one = cut_buck_global; if (narg == 6) cut_buck_one = force->numeric(FLERR,*(arg++)); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { buck_a_read[i][j] = buck_a_one; buck_c_read[i][j] = buck_c_one; buck_rho_read[i][j] = buck_rho_one; cut_buck_read[i][j] = cut_buck_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::init_style() { // require an atom style with charge defined if (!atom->q_flag && (ewald_order&(1<<1))) error->all(FLERR,"Pair style buck/long/coul/long requires atom attribute q"); // request regular or rRESPA neighbor lists if neighrequest_flag != 0 if (force->kspace->neighrequest_flag) { int irequest; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; if (respa == 0) irequest = neighbor->request(this,instance_me); else if (respa == 1) { irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } else { irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 2; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respamiddle = 1; irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } } else irequest = neighbor->request(this,instance_me); } cut_coulsq = cut_coul * cut_coul; // set rRESPA cutoffs if (strstr(update->integrate_style,"respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; else cut_respa = NULL; // ensure use of KSpace long-range solver, set two g_ewalds if (force->kspace == NULL) error->all(FLERR,"Pair style requires a KSpace style"); if (ewald_order&(1<<1)) g_ewald = force->kspace->g_ewald; if (ewald_order&(1<<6)) g_ewald_6 = force->kspace->g_ewald_6; // setup force tables - if (ncoultablebits) init_tables(cut_coul,cut_respa); - if (ndisptablebits) init_tables_disp(cut_buck_global); + if (ncoultablebits && (ewald_order&(1<<1))) init_tables(cut_coul,cut_respa); + if (ndisptablebits && (ewald_order&(1<<6))) init_tables_disp(cut_buck_global); } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use regular or rRESPA ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::init_list(int id, NeighList *ptr) { if (id == 0) list = ptr; else if (id == 1) listinner = ptr; else if (id == 2) listmiddle = ptr; else if (id == 3) listouter = ptr; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairBuckLongCoulLong::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); if (ewald_order&(1<<6)) cut_buck[i][j] = cut_buck_global; else cut_buck[i][j] = cut_buck_read[i][j]; buck_a[i][j] = buck_a_read[i][j]; buck_c[i][j] = buck_c_read[i][j]; buck_rho[i][j] = buck_rho_read[i][j]; double cut = MAX(cut_buck[i][j],cut_coul); cutsq[i][j] = cut*cut; cut_bucksq[i][j] = cut_buck[i][j] * cut_buck[i][j]; buck1[i][j] = buck_a[i][j]/buck_rho[i][j]; buck2[i][j] = 6.0*buck_c[i][j]; rhoinv[i][j] = 1.0/buck_rho[i][j]; // check interior rRESPA cutoff if (cut_respa && MIN(cut_buck[i][j],cut_coul) < cut_respa[3]) error->all(FLERR,"Pair cutoff < Respa interior cutoff"); if (offset_flag) { double rexp = exp(-cut_buck[i][j]/buck_rho[i][j]); offset[i][j] = buck_a[i][j]*rexp - buck_c[i][j]/pow(cut_buck[i][j],6.0); } else offset[i][j] = 0.0; cutsq[j][i] = cutsq[i][j]; cut_bucksq[j][i] = cut_bucksq[i][j]; buck_a[j][i] = buck_a[i][j]; buck_c[j][i] = buck_c[i][j]; rhoinv[j][i] = rhoinv[i][j]; buck1[j][i] = buck1[i][j]; buck2[j][i] = buck2[i][j]; offset[j][i] = offset[i][j]; return cut; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&buck_a_read[i][j],sizeof(double),1,fp); fwrite(&buck_rho_read[i][j],sizeof(double),1,fp); fwrite(&buck_c_read[i][j],sizeof(double),1,fp); fwrite(&cut_buck_read[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { fread(&buck_a_read[i][j],sizeof(double),1,fp); fread(&buck_rho_read[i][j],sizeof(double),1,fp); fread(&buck_c_read[i][j],sizeof(double),1,fp); fread(&cut_buck_read[i][j],sizeof(double),1,fp); } MPI_Bcast(&buck_a_read[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&buck_rho_read[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&buck_c_read[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_buck_read[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::write_restart_settings(FILE *fp) { fwrite(&cut_buck_global,sizeof(double),1,fp); fwrite(&cut_coul,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&ncoultablebits,sizeof(int),1,fp); fwrite(&tabinner,sizeof(double),1,fp); fwrite(&ewald_order,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&cut_buck_global,sizeof(double),1,fp); fread(&cut_coul,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); fread(&ncoultablebits,sizeof(int),1,fp); fread(&tabinner,sizeof(double),1,fp); fread(&ewald_order,sizeof(int),1,fp); } MPI_Bcast(&cut_buck_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); MPI_Bcast(&ewald_order,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g %g %g\n",i, buck_a_read[i][i],buck_rho_read[i][i],buck_c_read[i][i]); } /* ---------------------------------------------------------------------- proc 0 writes all pairs to data file ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) fprintf(fp,"%d %d %g %g %g\n",i,j, buck_a_read[i][j],buck_rho_read[i][j],buck_c_read[i][j]); } /* ---------------------------------------------------------------------- compute pair interactions ------------------------------------------------------------------------- */ void PairBuckLongCoulLong::compute(int eflag, int vflag) { double evdwl,ecoul,fpair; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x, *x0 = x[0]; double **f = atom->f, *f0 = f[0], *fi = f0; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; double qi = 0.0, qri = 0.0, *cutsqi, *cut_bucksqi, *buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti; double r, rsq, r2inv, force_coul, force_buck; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; vector xi, d; ineighn = (ineigh = list->ilist)+list->inum; for (; ineighfirstneigh[i])+list->numneigh[i]; for (; jneigh= cutsqi[typej = type[j]]) continue; r2inv = 1.0/rsq; r = sqrt(rsq); if (order1 && (rsq < cut_coulsq)) { // coulombic if (!ncoultablebits || rsq <= tabinnersq) { // series real space register double x = g_ewald*r; register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x); if (ni == 0) { s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s; if (eflag) ecoul = t; } else { // special case register double f = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f; if (eflag) ecoul = t-f; } } // table real space else { register union_int_float_t t; t.f = rsq; register const int k = (t.i & ncoulmask) >> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } else { // special case t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f); } } } else force_coul = ecoul = 0.0; if (rsq < cut_bucksqi[typej]) { // buckingham register double rn = r2inv*r2inv*r2inv, expr = exp(-r*rhoinvi[typej]); if (order6) { // long-range if (!ndisptablebits || rsq <= tabinnerdispsq) { register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*buckci[typej]; if (ni == 0) { force_buck = r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2; } else { // special case register double f = special_lj[ni], t = rn*(1.0-f); force_buck = f*r*expr*buck1i[typej]- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]; if (eflag) evdwl = f*expr*buckai[typej] - g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej]; } } else { //table real space register union_int_float_t disp_t; disp_t.f = rsq; register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits; register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k]; if (ni == 0) { force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]; if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]; } else { //speial case register double f = special_lj[ni], t = rn*(1.0-f); force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej]; if (eflag) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej]; } } } else { // cut if (ni == 0) { force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]; if (eflag) evdwl = expr*buckai[typej] - rn*buckci[typej]-offseti[typej]; } else { // special case register double f = special_lj[ni]; force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]); if (eflag) evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]); } } } else force_buck = evdwl = 0.0; fpair = (force_coul+force_buck)*r2inv; if (newton_pair || j < nlocal) { register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,ecoul,fpair,d[0],d[1],d[2]); } } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- */ void PairBuckLongCoulLong::compute_inner() { double r, rsq, r2inv, force_coul = 0.0, force_buck, fpair; int *type = atom->type; int nlocal = atom->nlocal; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi; vector xi, d; ineighn = (ineigh = listinner->ilist) + listinner->inum; for (; ineighfirstneigh[i])+listinner->numneigh[i]; for (; jneigh= cut_out_off_sq) continue; r2inv = 1.0/rsq; r = sqrt(rsq); if (order1 && (rsq < cut_coulsq)) // coulombic force_coul = ni == 0 ? qri*q[j]/r : qri*q[j]/r*special_coul[ni]; if (rsq < cut_bucksqi[typej = type[j]]) { // buckingham register double rn = r2inv*r2inv*r2inv, expr = exp(-r*rhoinvi[typej]); force_buck = ni == 0 ? (r*expr*buck1i[typej]-rn*buck2i[typej]) : (r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; } else force_buck = 0.0; fpair = (force_coul + force_buck) * r2inv; if (rsq > cut_out_on_sq) { // switching register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); } if (newton_pair || j < nlocal) { // force update register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } } } } /* ---------------------------------------------------------------------- */ void PairBuckLongCoulLong::compute_middle() { double r, rsq, r2inv, force_coul = 0.0, force_buck, fpair; int *type = atom->type; int nlocal = atom->nlocal; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; double cut_out_off = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); double qri, *cut_bucksqi, *buck1i, *buck2i, *rhoinvi; vector xi, d; ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum; for (; ineighfirstneigh[i])+listmiddle->numneigh[i]; for (; jneigh= cut_out_off_sq) continue; if (rsq <= cut_in_off_sq) continue; r2inv = 1.0/rsq; r = sqrt(rsq); if (order1 && (rsq < cut_coulsq)) // coulombic force_coul = ni == 0 ? qri*q[j]/r : qri*q[j]/r*special_coul[ni]; if (rsq < cut_bucksqi[typej = type[j]]) { // buckingham register double rn = r2inv*r2inv*r2inv, expr = exp(-r*rhoinvi[typej]); force_buck = ni == 0 ? (r*expr*buck1i[typej]-rn*buck2i[typej]) : (r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; } else force_buck = 0.0; fpair = (force_coul + force_buck) * r2inv; if (rsq < cut_in_on_sq) { // switching register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } if (rsq > cut_out_on_sq) { register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); } if (newton_pair || j < nlocal) { // force update register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } } } } /* ---------------------------------------------------------------------- */ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag) { double evdwl,ecoul,fpair,fvirial; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; double **x = atom->x, *x0 = x[0]; double **f = atom->f, *f0 = f[0], *fi = f0; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni, respa_flag; double qi = 0.0, qri = 0.0, *cutsqi, *cut_bucksqi, *buck1i, *buck2i, *buckai, *buckci, *rhoinvi, *offseti; double r, rsq, r2inv, force_coul, force_buck; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; double respa_buck = 0.0, respa_coul = 0.0, frespa = 0.0; vector xi, d; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; ineighn = (ineigh = listouter->ilist)+listouter->inum; for (; ineighfirstneigh[i])+listouter->numneigh[i]; for (; jneigh= cutsqi[typej = type[j]]) continue; r2inv = 1.0/rsq; r = sqrt(rsq); frespa = 1.0; //check whether and how to compute respa corrections respa_coul = 0.0; respa_buck = 0.0; respa_flag = rsq < cut_in_on_sq ? 1 : 0; if (respa_flag && (rsq > cut_in_off_sq)) { register double rsw = (r-cut_in_off)/cut_in_diff; frespa = 1-rsw*rsw*(3.0-2.0*rsw); } if (order1 && (rsq < cut_coulsq)) { // coulombic if (!ncoultablebits || rsq <= tabinnersq) { // series real space register double s = qri*q[j]; if (respa_flag) // correct for respa respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni]; register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x); if (ni == 0) { s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul; if (eflag) ecoul = t; } else { // correct for special register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul; if (eflag) ecoul = t-ri; } } // table real space else { if (respa_flag) { register double s = qri*q[j]; respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni]; } register union_int_float_t t; t.f = rsq; register const int k = (t.i & ncoulmask) >> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } else { // correct for special t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); if (eflag) { t.f = (1.0-special_coul[ni])*(ptable[k]+f*dptable[k]); ecoul = qiqj*(etable[k]+f*detable[k]-t.f); } } } } else force_coul = respa_coul = ecoul = 0.0; if (rsq < cut_bucksqi[typej]) { // buckingham register double rn = r2inv*r2inv*r2inv, expr = exp(-r*rhoinvi[typej]); if (respa_flag) respa_buck = ni == 0 ? // correct for respa frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) : frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; if (order6) { // long-range form if (!ndisptablebits || rsq <= tabinnerdispsq) { register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*buckci[typej]; if (ni == 0) { force_buck = r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq-respa_buck; if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2; } else { // correct for special register double f = special_lj[ni], t = rn*(1.0-f); force_buck = f*r*expr*buck1i[typej]- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck; if (eflag) evdwl = f*expr*buckai[typej] - g6*((a2+1.0)*a2+0.5)*x2+t*buckci[typej]; } } else { // table real space register union_int_float_t disp_t; disp_t.f = rsq; register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits; register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k]; register double rn = r2inv*r2inv*r2inv; if (ni == 0) { force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck; if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]; } else { //special case register double f = special_lj[ni], t = rn*(1.0-f); force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck; if (eflag) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej]; } } } else { // cut form if (ni == 0) { force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]-respa_buck; if (eflag) evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej]; } else { // correct for special register double f = special_lj[ni]; force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck; if (eflag) evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]); } } } else force_buck = respa_buck = evdwl = 0.0; fpair = (force_coul+force_buck)*r2inv; if (newton_pair || j < nlocal) { register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } if (evflag) { fvirial = (force_coul + force_buck + respa_coul + respa_buck)*r2inv; ev_tally(i,j,nlocal,newton_pair, evdwl,ecoul,fvirial,d[0],d[1],d[2]); } } } } /* ---------------------------------------------------------------------- */ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_buck, double &fforce) { double f, r, r2inv, r6inv, force_coul, force_buck; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q; r = sqrt(rsq); r2inv = 1.0/rsq; double eng = 0.0; if ((ewald_order&2) && (rsq < cut_coulsq)) { // coulombic if (!ncoultablebits || rsq <= tabinnersq) { // series real space register double x = g_ewald*r; register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x); f = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f; eng += t-f; } else { // table real space register union_int_float_t t; t.f = rsq; register const int k = (t.i & ncoulmask) >> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j]; t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]); force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); eng += qiqj*(etable[k]+f*detable[k]-t.f); } } else force_coul = 0.0; if (rsq < cut_bucksq[itype][jtype]) { // buckingham register double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]); r6inv = r2inv*r2inv*r2inv; if (ewald_order&64) { // long-range register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_buck); x2 = a2*exp(-x2)*buck_c[itype][jtype]; force_buck = buck1[itype][jtype]*r*expr- g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*buck2[itype][jtype]; eng += buck_a[itype][jtype]*expr- g6*((a2+1.0)*a2+0.5)*x2+t*buck_c[itype][jtype]; } else { // cut force_buck = buck1[itype][jtype]*r*expr-factor_buck*buck_c[itype][jtype]*r6inv; eng += buck_a[itype][jtype]*expr- factor_buck*(buck_c[itype][jtype]*r6inv-offset[itype][jtype]); } } else force_buck = 0.0; fforce = (force_coul+force_buck)*r2inv; return eng; } diff --git a/src/KSPACE/pair_lj_long_coul_long.cpp b/src/KSPACE/pair_lj_long_coul_long.cpp index 3675ea76d..8197d16f7 100644 --- a/src/KSPACE/pair_lj_long_coul_long.cpp +++ b/src/KSPACE/pair_lj_long_coul_long.cpp @@ -1,1041 +1,1042 @@ /* ---------------------------------------------------------------------- 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: Pieter J. in 't Veld (SNL) Tabulation for long-range dispersion added by Wayne Mitchell (Loyola University New Orleans) ------------------------------------------------------------------------- */ #include #include #include #include #include "math_vector.h" #include "pair_lj_long_coul_long.h" #include "atom.h" #include "comm.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "force.h" #include "kspace.h" #include "update.h" #include "integrate.h" #include "respa.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; #define EWALD_F 1.12837917 #define EWALD_P 0.3275911 #define A1 0.254829592 #define A2 -0.284496736 #define A3 1.421413741 #define A4 -1.453152027 #define A5 1.061405429 /* ---------------------------------------------------------------------- */ PairLJLongCoulLong::PairLJLongCoulLong(LAMMPS *lmp) : Pair(lmp) { dispersionflag = ewaldflag = pppmflag = 1; respa_enable = 1; writedata = 1; ftable = NULL; fdisptable = NULL; qdist = 0.0; } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairLJLongCoulLong::options(char **arg, int order) { const char *option[] = {"long", "cut", "off", NULL}; int i; if (!*arg) error->all(FLERR,"Illegal pair_style lj/long/coul/long command"); for (i=0; option[i]&&strcmp(arg[0], option[i]); ++i); switch (i) { default: error->all(FLERR,"Illegal pair_style lj/long/coul/long command"); case 0: ewald_order |= 1<all(FLERR,"Illegal pair_style command"); ewald_off = 0; ewald_order = 0; options(arg, 6); options(++arg, 1); if (!comm->me && ewald_order == ((1<<1) | (1<<6))) error->warning(FLERR,"Using largest cutoff for lj/long/coul/long"); if (!*(++arg)) error->all(FLERR,"Cutoffs missing in pair_style lj/long/coul/long"); if (!((ewald_order^ewald_off) & (1<<1))) error->all(FLERR, "Coulomb cut not supported in pair_style lj/long/coul/long"); cut_lj_global = force->numeric(FLERR,*(arg++)); if (narg == 4 && ((ewald_order & 0x42) == 0x42)) error->all(FLERR,"Only one cutoff allowed when requesting all long"); if (narg == 4) cut_coul = force->numeric(FLERR,*arg); else cut_coul = cut_lj_global; if (allocated) { int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i+1; j <= atom->ntypes; j++) if (setflag[i][j]) cut_lj[i][j] = cut_lj_global; } } /* ---------------------------------------------------------------------- free all arrays ------------------------------------------------------------------------- */ PairLJLongCoulLong::~PairLJLongCoulLong() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(cut_lj_read); memory->destroy(cut_lj); memory->destroy(cut_ljsq); memory->destroy(epsilon_read); memory->destroy(epsilon); memory->destroy(sigma_read); memory->destroy(sigma); memory->destroy(lj1); memory->destroy(lj2); memory->destroy(lj3); memory->destroy(lj4); memory->destroy(offset); } if (ftable) free_tables(); + if (fdisptable) free_disp_tables(); } /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ void PairLJLongCoulLong::allocate() { allocated = 1; int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); for (int i = 1; i <= n; i++) for (int j = i; j <= n; j++) setflag[i][j] = 0; memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cut_lj_read,n+1,n+1,"pair:cut_lj_read"); memory->create(cut_lj,n+1,n+1,"pair:cut_lj"); memory->create(cut_ljsq,n+1,n+1,"pair:cut_ljsq"); memory->create(epsilon_read,n+1,n+1,"pair:epsilon_read"); memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(sigma_read,n+1,n+1,"pair:sigma_read"); memory->create(sigma,n+1,n+1,"pair:sigma"); memory->create(lj1,n+1,n+1,"pair:lj1"); memory->create(lj2,n+1,n+1,"pair:lj2"); memory->create(lj3,n+1,n+1,"pair:lj3"); memory->create(lj4,n+1,n+1,"pair:lj4"); memory->create(offset,n+1,n+1,"pair:offset"); } /* ---------------------------------------------------------------------- extract protected data from object ------------------------------------------------------------------------- */ void *PairLJLongCoulLong::extract(const char *id, int &dim) { const char *ids[] = { "B", "sigma", "epsilon", "ewald_order", "ewald_cut", "ewald_mix", "cut_coul", "cut_LJ", NULL}; void *ptrs[] = { lj4, sigma, epsilon, &ewald_order, &cut_coul, &mix_flag, &cut_coul, &cut_lj_global, NULL}; int i; for (i=0; ids[i]&&strcmp(ids[i], id); ++i); if (i <= 2) dim = 2; else dim = 0; return ptrs[i]; } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairLJLongCoulLong::coeff(int narg, char **arg) { if (narg < 4 || narg > 5) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); int ilo,ihi,jlo,jhi; force->bounds(arg[0],atom->ntypes,ilo,ihi); force->bounds(arg[1],atom->ntypes,jlo,jhi); double epsilon_one = force->numeric(FLERR,arg[2]); double sigma_one = force->numeric(FLERR,arg[3]); double cut_lj_one = cut_lj_global; if (narg == 5) cut_lj_one = force->numeric(FLERR,arg[4]); int count = 0; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { epsilon_read[i][j] = epsilon_one; sigma_read[i][j] = sigma_one; cut_lj_read[i][j] = cut_lj_one; setflag[i][j] = 1; count++; } } if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairLJLongCoulLong::init_style() { // require an atom style with charge defined if (!atom->q_flag && (ewald_order&(1<<1))) error->all(FLERR, "Invoking coulombic in pair style lj/coul requires atom attribute q"); // request regular or rRESPA neighbor lists if neighrequest_flag != 0 if (force->kspace->neighrequest_flag) { int irequest; if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) { int respa = 0; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; if (respa == 0) irequest = neighbor->request(this,instance_me); else if (respa == 1) { irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } else { irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 1; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respainner = 1; irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 2; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respamiddle = 1; irequest = neighbor->request(this,instance_me); neighbor->requests[irequest]->id = 3; neighbor->requests[irequest]->half = 0; neighbor->requests[irequest]->respaouter = 1; } } else irequest = neighbor->request(this,instance_me); } cut_coulsq = cut_coul * cut_coul; // set rRESPA cutoffs if (strstr(update->integrate_style,"respa") && ((Respa *) update->integrate)->level_inner >= 0) cut_respa = ((Respa *) update->integrate)->cutoff; else cut_respa = NULL; // ensure use of KSpace long-range solver, set g_ewald if (force->kspace == NULL) error->all(FLERR,"Pair style requires a KSpace style"); if (force->kspace) g_ewald = force->kspace->g_ewald; if (force->kspace) g_ewald_6 = force->kspace->g_ewald_6; // setup force tables - if (ncoultablebits) init_tables(cut_coul,cut_respa); - if (ndisptablebits) init_tables_disp(cut_lj_global); + if (ncoultablebits && (ewald_order&(1<<1))) init_tables(cut_coul,cut_respa); + if (ndisptablebits && (ewald_order&(1<<6))) init_tables_disp(cut_lj_global); } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use regular or rRESPA ------------------------------------------------------------------------- */ void PairLJLongCoulLong::init_list(int id, NeighList *ptr) { if (id == 0) list = ptr; else if (id == 1) listinner = ptr; else if (id == 2) listmiddle = ptr; else if (id == 3) listouter = ptr; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairLJLongCoulLong::init_one(int i, int j) { if (setflag[i][j] == 0) { epsilon[i][j] = mix_energy(epsilon_read[i][i],epsilon_read[j][j], sigma_read[i][i],sigma_read[j][j]); sigma[i][j] = mix_distance(sigma_read[i][i],sigma_read[j][j]); if (ewald_order&(1<<6)) cut_lj[i][j] = cut_lj_global; else cut_lj[i][j] = mix_distance(cut_lj_read[i][i],cut_lj_read[j][j]); } else { sigma[i][j] = sigma_read[i][j]; epsilon[i][j] = epsilon_read[i][j]; cut_lj[i][j] = cut_lj_read[i][j]; } double cut = MAX(cut_lj[i][j], cut_coul + 2.0*qdist); cutsq[i][j] = cut*cut; cut_ljsq[i][j] = cut_lj[i][j] * cut_lj[i][j]; lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); // check interior rRESPA cutoff if (cut_respa && MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) error->all(FLERR,"Pair cutoff < Respa interior cutoff"); if (offset_flag) { double ratio = sigma[i][j] / cut_lj[i][j]; offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); } else offset[i][j] = 0.0; cutsq[j][i] = cutsq[i][j]; cut_ljsq[j][i] = cut_ljsq[i][j]; lj1[j][i] = lj1[i][j]; lj2[j][i] = lj2[i][j]; lj3[j][i] = lj3[i][j]; lj4[j][i] = lj4[i][j]; offset[j][i] = offset[i][j]; return cut; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJLongCoulLong::write_restart(FILE *fp) { write_restart_settings(fp); int i,j; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { fwrite(&setflag[i][j],sizeof(int),1,fp); if (setflag[i][j]) { fwrite(&epsilon_read[i][j],sizeof(double),1,fp); fwrite(&sigma_read[i][j],sizeof(double),1,fp); fwrite(&cut_lj_read[i][j],sizeof(double),1,fp); } } } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJLongCoulLong::read_restart(FILE *fp) { read_restart_settings(fp); allocate(); int i,j; int me = comm->me; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); if (setflag[i][j]) { if (me == 0) { fread(&epsilon_read[i][j],sizeof(double),1,fp); fread(&sigma_read[i][j],sizeof(double),1,fp); fread(&cut_lj_read[i][j],sizeof(double),1,fp); } MPI_Bcast(&epsilon_read[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma_read[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_lj_read[i][j],1,MPI_DOUBLE,0,world); } } } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairLJLongCoulLong::write_restart_settings(FILE *fp) { fwrite(&cut_lj_global,sizeof(double),1,fp); fwrite(&cut_coul,sizeof(double),1,fp); fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&ncoultablebits,sizeof(int),1,fp); fwrite(&tabinner,sizeof(double),1,fp); fwrite(&ewald_order,sizeof(int),1,fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairLJLongCoulLong::read_restart_settings(FILE *fp) { if (comm->me == 0) { fread(&cut_lj_global,sizeof(double),1,fp); fread(&cut_coul,sizeof(double),1,fp); fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); fread(&ncoultablebits,sizeof(int),1,fp); fread(&tabinner,sizeof(double),1,fp); fread(&ewald_order,sizeof(int),1,fp); } MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); MPI_Bcast(&ewald_order,1,MPI_INT,0,world); } /* ---------------------------------------------------------------------- proc 0 writes to data file ------------------------------------------------------------------------- */ void PairLJLongCoulLong::write_data(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g %g\n",i,epsilon_read[i][i],sigma_read[i][i]); } /* ---------------------------------------------------------------------- proc 0 writes all pairs to data file ------------------------------------------------------------------------- */ void PairLJLongCoulLong::write_data_all(FILE *fp) { for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) fprintf(fp,"%d %d %g %g %g\n",i,j, epsilon_read[i][j],sigma_read[i][j],cut_lj_read[i][j]); } /* ---------------------------------------------------------------------- compute pair interactions ------------------------------------------------------------------------- */ void PairLJLongCoulLong::compute(int eflag, int vflag) { double evdwl,ecoul,fpair; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x, *x0 = x[0]; double **f = atom->f, *f0 = f[0], *fi = f0; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; double qi = 0.0, qri = 0.0; double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti; double rsq, r2inv, force_coul, force_lj; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; vector xi, d; ineighn = (ineigh = list->ilist)+list->inum; for (; ineighfirstneigh[i])+list->numneigh[i]; for (; jneigh= cutsqi[typej = type[j]]) continue; r2inv = 1.0/rsq; if (order1 && (rsq < cut_coulsq)) { // coulombic if (!ncoultablebits || rsq <= tabinnersq) { // series real space register double r = sqrt(rsq), x = g_ewald*r; register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x); if (ni == 0) { s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s; if (eflag) ecoul = t; } else { // special case r = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r; if (eflag) ecoul = t-r; } } // table real space else { register union_int_float_t t; t.f = rsq; register const int k = (t.i & ncoulmask)>>ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } else { // special case t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f); } } } else force_coul = ecoul = 0.0; if (rsq < cut_ljsqi[typej]) { // lj if (order6) { // long-range lj if(!ndisptablebits || rsq <= tabinnerdispsq) { // series real space register double rn = r2inv*r2inv*r2inv; register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*lj4i[typej]; if (ni == 0) { force_lj = (rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2; } else { // special case register double f = special_lj[ni], t = rn*(1.0-f); force_lj = f*(rn *= rn)*lj1i[typej]- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]; if (eflag) evdwl = f*rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[typej]; } } else { // table real space register union_int_float_t disp_t; disp_t.f = rsq; register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits; register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k]; register double rn = r2inv*r2inv*r2inv; if (ni == 0) { force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]; if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]; } else { // special case register double f = special_lj[ni], t = rn*(1.0-f); force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]; if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej]; } } } else { // cut lj register double rn = r2inv*r2inv*r2inv; if (ni == 0) { force_lj = rn*(rn*lj1i[typej]-lj2i[typej]); if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]; } else { // special case register double f = special_lj[ni]; force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]); if (eflag) evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]); } } } else force_lj = evdwl = 0.0; fpair = (force_coul+force_lj)*r2inv; if (newton_pair || j < nlocal) { register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } if (evflag) ev_tally(i,j,nlocal,newton_pair, evdwl,ecoul,fpair,d[0],d[1],d[2]); } } if (vflag_fdotr) virial_fdotr_compute(); } /* ---------------------------------------------------------------------- */ void PairLJLongCoulLong::compute_inner() { double rsq, r2inv, force_coul = 0.0, force_lj, fpair; int *type = atom->type; int nlocal = atom->nlocal; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; double cut_out_on = cut_respa[0]; double cut_out_off = cut_respa[1]; double cut_out_diff = cut_out_off - cut_out_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); double qri, *cut_ljsqi, *lj1i, *lj2i; vector xi, d; ineighn = (ineigh = listinner->ilist)+listinner->inum; for (; ineighfirstneigh[i])+listinner->numneigh[i]; for (; jneigh= cut_out_off_sq) continue; r2inv = 1.0/rsq; if (order1 && (rsq < cut_coulsq)) { // coulombic qri = qqrd2e*q[i]; force_coul = ni == 0 ? qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni]; } if (rsq < cut_ljsqi[typej = type[j]]) { // lennard-jones register double rn = r2inv*r2inv*r2inv; force_lj = ni == 0 ? rn*(rn*lj1i[typej]-lj2i[typej]) : rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; } else force_lj = 0.0; fpair = (force_coul + force_lj) * r2inv; if (rsq > cut_out_on_sq) { // switching register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); } if (newton_pair || j < nlocal) { // force update register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } } } } /* ---------------------------------------------------------------------- */ void PairLJLongCoulLong::compute_middle() { double rsq, r2inv, force_coul = 0.0, force_lj, fpair; int *type = atom->type; int nlocal = atom->nlocal; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; double cut_in_off = cut_respa[0]; double cut_in_on = cut_respa[1]; double cut_out_on = cut_respa[2]; double cut_out_off = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_out_diff = cut_out_off - cut_out_on; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_off_sq = cut_out_off*cut_out_off; int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni; int i, j, order1 = (ewald_order|(ewald_off^-1))&(1<<1); double qri, *cut_ljsqi, *lj1i, *lj2i; vector xi, d; ineighn = (ineigh = listmiddle->ilist)+listmiddle->inum; for (; ineighfirstneigh[i])+listmiddle->numneigh[i]; for (; jneigh= cut_out_off_sq) continue; if (rsq <= cut_in_off_sq) continue; r2inv = 1.0/rsq; if (order1 && (rsq < cut_coulsq)) // coulombic force_coul = ni == 0 ? qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni]; if (rsq < cut_ljsqi[typej = type[j]]) { // lennard-jones register double rn = r2inv*r2inv*r2inv; force_lj = ni == 0 ? rn*(rn*lj1i[typej]-lj2i[typej]) : rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; } else force_lj = 0.0; fpair = (force_coul + force_lj) * r2inv; if (rsq < cut_in_on_sq) { // switching register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; fpair *= rsw*rsw*(3.0 - 2.0*rsw); } if (rsq > cut_out_on_sq) { register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0); } if (newton_pair || j < nlocal) { // force update register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } } } } /* ---------------------------------------------------------------------- */ void PairLJLongCoulLong::compute_outer(int eflag, int vflag) { double evdwl,ecoul,fvirial,fpair; evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; double **x = atom->x, *x0 = x[0]; double **f = atom->f, *f0 = f[0], *fi = f0; double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; int i, j, order1 = ewald_order&(1<<1), order6 = ewald_order&(1<<6); int *ineigh, *ineighn, *jneigh, *jneighn, typei, typej, ni, respa_flag; double qi = 0.0, qri = 0.0; double *cutsqi, *cut_ljsqi, *lj1i, *lj2i, *lj3i, *lj4i, *offseti; double rsq, r2inv, force_coul, force_lj; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; double respa_lj = 0.0, respa_coul = 0.0, frespa = 0.0; vector xi, d; double cut_in_off = cut_respa[2]; double cut_in_on = cut_respa[3]; double cut_in_diff = cut_in_on - cut_in_off; double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_on_sq = cut_in_on*cut_in_on; ineighn = (ineigh = listouter->ilist)+listouter->inum; for (; ineighfirstneigh[i])+listouter->numneigh[i]; for (; jneigh= cutsqi[typej = type[j]]) continue; r2inv = 1.0/rsq; frespa = 1.0; // check whether and how to compute respa corrections respa_coul = 0; respa_lj = 0; respa_flag = rsq < cut_in_on_sq ? 1 : 0; if (respa_flag && (rsq > cut_in_off_sq)) { register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff; frespa = 1-rsw*rsw*(3.0-2.0*rsw); } if (order1 && (rsq < cut_coulsq)) { // coulombic if (!ncoultablebits || rsq <= tabinnersq) { // series real space register double r = sqrt(rsq), s = qri*q[j]; if (respa_flag) // correct for respa respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni]; register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x); if (ni == 0) { s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul; if (eflag) ecoul = t; } else { // correct for special r = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r-respa_coul; if (eflag) ecoul = t-r; } } // table real space else { if (respa_flag) { register double r = sqrt(rsq), s = qri*q[j]; respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni]; } register union_int_float_t t; t.f = rsq; register const int k = (t.i & ncoulmask) >> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } else { // correct for special t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]); force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); if (eflag) { t.f = (1.0-special_coul[ni])*(ptable[k]+f*dptable[k]); ecoul = qiqj*(etable[k]+f*detable[k]-t.f); } } } } else force_coul = respa_coul = ecoul = 0.0; if (rsq < cut_ljsqi[typej]) { // lennard-jones register double rn = r2inv*r2inv*r2inv; if (respa_flag) respa_lj = ni == 0 ? // correct for respa frespa*rn*(rn*lj1i[typej]-lj2i[typej]) : frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; if (order6) { // long-range form if (!ndisptablebits || rsq <= tabinnerdispsq) { register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*lj4i[typej]; if (ni == 0) { force_lj = (rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq-respa_lj; if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2; } else { // correct for special register double f = special_lj[ni], t = rn*(1.0-f); force_lj = f*(rn *= rn)*lj1i[typej]- g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj; if (eflag) evdwl = f*rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2+t*lj4i[typej]; } } else { // table real space register union_int_float_t disp_t; disp_t.f = rsq; register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits; register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k]; register double rn = r2inv*r2inv*r2inv; if (ni == 0) { force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj; if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]; } else { // special case register double f = special_lj[ni], t = rn*(1.0-f); force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj; if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej]; } } } else { // cut form if (ni == 0) { force_lj = rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj; if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]; } else { // correct for special register double f = special_lj[ni]; force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj; if (eflag) evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]); } } } else force_lj = respa_lj = evdwl = 0.0; fpair = (force_coul+force_lj)*r2inv; if (newton_pair || j < nlocal) { register double *fj = f0+(j+(j<<1)), f; fi[0] += f = d[0]*fpair; fj[0] -= f; fi[1] += f = d[1]*fpair; fj[1] -= f; fi[2] += f = d[2]*fpair; fj[2] -= f; } else { fi[0] += d[0]*fpair; fi[1] += d[1]*fpair; fi[2] += d[2]*fpair; } if (evflag) { fvirial = (force_coul + force_lj + respa_coul + respa_lj)*r2inv; ev_tally(i,j,nlocal,newton_pair, evdwl,ecoul,fvirial,d[0],d[1],d[2]); } } } } /* ---------------------------------------------------------------------- */ double PairLJLongCoulLong::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { double r2inv, r6inv, force_coul, force_lj; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2, *q = atom->q; double eng = 0.0; r2inv = 1.0/rsq; if ((ewald_order&2) && (rsq < cut_coulsq)) { // coulombic if (!ncoultablebits || rsq <= tabinnersq) { // series real space register double r = sqrt(rsq), x = g_ewald*r; register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x); r = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x); force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r; eng += t-r; } else { // table real space register union_int_float_t t; t.f = rsq; register const int k = (t.i & ncoulmask) >> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j]; t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]); force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f); eng += qiqj*(etable[k]+f*detable[k]-t.f); } } else force_coul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { // lennard-jones r6inv = r2inv*r2inv*r2inv; if (ewald_order&64) { // long-range register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj); x2 = a2*exp(-x2)*lj4[itype][jtype]; force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]- g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype]; eng += factor_lj*r6inv*lj3[itype][jtype]- g6*((a2+1.0)*a2+0.5)*x2+t*lj4[itype][jtype]; } else { // cut force_lj = factor_lj*r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]); eng += factor_lj*(r6inv*(r6inv*lj3[itype][jtype]- lj4[itype][jtype])-offset[itype][jtype]); } } else force_lj = 0.0; fforce = (force_coul+force_lj)*r2inv; return eng; } diff --git a/src/pair.cpp b/src/pair.cpp index 651cabed6..5d73a592e 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -1,1717 +1,1719 @@ /* ---------------------------------------------------------------------- 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: Paul Crozier (SNL) ------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include "pair.h" #include "atom.h" #include "neighbor.h" #include "neigh_list.h" #include "domain.h" #include "comm.h" #include "force.h" #include "kspace.h" #include "update.h" #include "modify.h" #include "compute.h" #include "suffix.h" #include "atom_masks.h" #include "memory.h" #include "math_const.h" #include "error.h" using namespace LAMMPS_NS; using namespace MathConst; enum{NONE,RLINEAR,RSQ,BMP}; // allocate space for static class instance variable and initialize it int Pair::instance_total = 0; /* ---------------------------------------------------------------------- */ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) { instance_me = instance_total++; eng_vdwl = eng_coul = 0.0; comm_forward = comm_reverse = comm_reverse_off = 0; single_enable = 1; restartinfo = 1; respa_enable = 0; one_coeff = 0; no_virial_fdotr_compute = 0; writedata = 0; ghostneigh = 0; nextra = 0; pvector = NULL; single_extra = 0; svector = NULL; ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag = dipoleflag = 0; reinitflag = 1; // pair_modify settingsx compute_flag = 1; manybody_flag = 0; offset_flag = 0; mix_flag = GEOMETRIC; tail_flag = 0; etail = ptail = etail_ij = ptail_ij = 0.0; ncoultablebits = 12; ndisptablebits = 12; tabinner = sqrt(2.0); tabinner_disp = sqrt(2.0); + ftable = NULL; + fdisptable = NULL; allocated = 0; suffix_flag = Suffix::NONE; maxeatom = maxvatom = 0; eatom = NULL; vatom = NULL; num_tally_compute = 0; list_tally_compute = NULL; // KOKKOS per-fix data masks execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; copymode = 0; } /* ---------------------------------------------------------------------- */ Pair::~Pair() { num_tally_compute = 0; memory->sfree((void *) list_tally_compute); list_tally_compute = NULL; if (copymode) return; memory->destroy(eatom); memory->destroy(vatom); } /* ---------------------------------------------------------------------- modify parameters of the pair style pair_hybrid has its own version of this routine to apply modifications to each of its sub-styles ------------------------------------------------------------------------- */ void Pair::modify_params(int narg, char **arg) { if (narg == 0) error->all(FLERR,"Illegal pair_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"mix") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"geometric") == 0) mix_flag = GEOMETRIC; else if (strcmp(arg[iarg+1],"arithmetic") == 0) mix_flag = ARITHMETIC; else if (strcmp(arg[iarg+1],"sixthpower") == 0) mix_flag = SIXTHPOWER; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"shift") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) offset_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) offset_flag = 0; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"table") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); ncoultablebits = force->inumeric(FLERR,arg[iarg+1]); if (ncoultablebits > sizeof(float)*CHAR_BIT) error->all(FLERR,"Too many total bits for bitmapped lookup table"); iarg += 2; } else if (strcmp(arg[iarg],"table/disp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); ndisptablebits = force->inumeric(FLERR,arg[iarg+1]); if (ndisptablebits > sizeof(float)*CHAR_BIT) error->all(FLERR,"Too many total bits for bitmapped lookup table"); iarg += 2; } else if (strcmp(arg[iarg],"tabinner") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); tabinner = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"tabinner/disp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); tabinner_disp = force->numeric(FLERR,arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"tail") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) tail_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) tail_flag = 0; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else if (strcmp(arg[iarg],"compute") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal pair_modify command"); if (strcmp(arg[iarg+1],"yes") == 0) compute_flag = 1; else if (strcmp(arg[iarg+1],"no") == 0) compute_flag = 0; else error->all(FLERR,"Illegal pair_modify command"); iarg += 2; } else error->all(FLERR,"Illegal pair_modify command"); } } /* ---------------------------------------------------------------------- */ void Pair::init() { int i,j; if (offset_flag && tail_flag) error->all(FLERR,"Cannot have both pair_modify shift and tail set to yes"); if (tail_flag && domain->dimension == 2) error->all(FLERR,"Cannot use pair tail corrections with 2d simulations"); if (tail_flag && domain->nonperiodic && comm->me == 0) error->warning(FLERR,"Using pair tail corrections with nonperiodic system"); if (!compute_flag && tail_flag) error->warning(FLERR,"Using pair tail corrections with " "pair_modify compute no"); if (!compute_flag && offset_flag) error->warning(FLERR,"Using pair potential shift with " "pair_modify compute no"); // for manybody potentials // check if bonded exclusions could invalidate the neighbor list if (manybody_flag && atom->molecular) { int flag = 0; if (atom->nbonds > 0 && force->special_lj[1] == 0.0 && force->special_coul[1] == 0.0) flag = 1; if (atom->nangles > 0 && force->special_lj[2] == 0.0 && force->special_coul[2] == 0.0) flag = 1; if (atom->ndihedrals > 0 && force->special_lj[3] == 0.0 && force->special_coul[3] == 0.0) flag = 1; if (flag && comm->me == 0) error->warning(FLERR,"Using a manybody potential with " "bonds/angles/dihedrals and special_bond exclusions"); } // I,I coeffs must be set // init_one() will check if I,J is set explicitly or inferred by mixing if (!allocated) error->all(FLERR,"All pair coeffs are not set"); for (i = 1; i <= atom->ntypes; i++) if (setflag[i][i] == 0) error->all(FLERR,"All pair coeffs are not set"); // style-specific initialization init_style(); // call init_one() for each I,J // set cutsq for each I,J, used to neighbor // cutforce = max of all I,J cutoffs cutforce = 0.0; etail = ptail = 0.0; double cut; for (i = 1; i <= atom->ntypes; i++) for (j = i; j <= atom->ntypes; j++) { cut = init_one(i,j); cutsq[i][j] = cutsq[j][i] = cut*cut; cutforce = MAX(cutforce,cut); if (tail_flag) { etail += etail_ij; ptail += ptail_ij; if (i != j) { etail += etail_ij; ptail += ptail_ij; } } } } /* ---------------------------------------------------------------------- reset all type-based params by invoking init_one() for each I,J called by fix adapt after it changes one or more params ------------------------------------------------------------------------- */ void Pair::reinit() { // generalize this error message if reinit() is used by more than fix adapt if (!reinitflag) error->all(FLERR,"Fix adapt interface to this pair style not supported"); etail = ptail = 0.0; for (int i = 1; i <= atom->ntypes; i++) for (int j = i; j <= atom->ntypes; j++) { init_one(i,j); if (tail_flag) { etail += etail_ij; ptail += ptail_ij; if (i != j) { etail += etail_ij; ptail += ptail_ij; } } } } /* ---------------------------------------------------------------------- init specific to a pair style specific pair style can override this function if needs its own error checks if needs another kind of neighbor list request default neighbor list = half list ------------------------------------------------------------------------- */ void Pair::init_style() { neighbor->request(this,instance_me); } /* ---------------------------------------------------------------------- neighbor callback to inform pair style of neighbor list to use specific pair style can override this function ------------------------------------------------------------------------- */ void Pair::init_list(int which, NeighList *ptr) { list = ptr; } /* ---------------------------------------------------------------------- setup Coulomb force tables used in compute routines ------------------------------------------------------------------------- */ void Pair::init_tables(double cut_coul, double *cut_respa) { int masklo,maskhi; double r,grij,expm2,derfc,egamma,fgamma,rsw; double qqrd2e = force->qqrd2e; if (force->kspace == NULL) error->all(FLERR,"Pair style requires a KSpace style"); double g_ewald = force->kspace->g_ewald; double cut_coulsq = cut_coul * cut_coul; tabinnersq = tabinner*tabinner; init_bitmap(tabinner,cut_coul,ncoultablebits, masklo,maskhi,ncoulmask,ncoulshiftbits); int ntable = 1; for (int i = 0; i < ncoultablebits; i++) ntable *= 2; // linear lookup tables of length N = 2^ncoultablebits // stored value = value at lower edge of bin // d values = delta from lower edge to upper edge of bin if (ftable) free_tables(); memory->create(rtable,ntable,"pair:rtable"); memory->create(ftable,ntable,"pair:ftable"); memory->create(ctable,ntable,"pair:ctable"); memory->create(etable,ntable,"pair:etable"); memory->create(drtable,ntable,"pair:drtable"); memory->create(dftable,ntable,"pair:dftable"); memory->create(dctable,ntable,"pair:dctable"); memory->create(detable,ntable,"pair:detable"); if (cut_respa == NULL) { vtable = ptable = dvtable = dptable = NULL; } else { memory->create(vtable,ntable,"pair:vtable"); memory->create(ptable,ntable,"pair:ptable"); memory->create(dvtable,ntable,"pair:dvtable"); memory->create(dptable,ntable,"pair:dptable"); } union_int_float_t rsq_lookup; union_int_float_t minrsq_lookup; int itablemin; minrsq_lookup.i = 0 << ncoulshiftbits; minrsq_lookup.i |= maskhi; for (int i = 0; i < ntable; i++) { rsq_lookup.i = i << ncoulshiftbits; rsq_lookup.i |= masklo; if (rsq_lookup.f < tabinnersq) { rsq_lookup.i = i << ncoulshiftbits; rsq_lookup.i |= maskhi; } r = sqrtf(rsq_lookup.f); if (msmflag) { egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)* force->kspace->dgamma(r/cut_coul); } else { grij = g_ewald * r; expm2 = exp(-grij*grij); derfc = erfc(grij); } if (cut_respa == NULL) { rtable[i] = rsq_lookup.f; ctable[i] = qqrd2e/r; if (msmflag) { ftable[i] = qqrd2e/r * fgamma; etable[i] = qqrd2e/r * egamma; } else { ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); etable[i] = qqrd2e/r * derfc; } } else { rtable[i] = rsq_lookup.f; ctable[i] = 0.0; ptable[i] = qqrd2e/r; if (msmflag) { ftable[i] = qqrd2e/r * (fgamma - 1.0); etable[i] = qqrd2e/r * egamma; vtable[i] = qqrd2e/r * fgamma; } else { ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0); etable[i] = qqrd2e/r * derfc; vtable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); } if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); } else { if (msmflag) ftable[i] = qqrd2e/r * fgamma; else ftable[i] = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); ctable[i] = qqrd2e/r; } } } minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); } tabinnersq = minrsq_lookup.f; int ntablem1 = ntable - 1; for (int i = 0; i < ntablem1; i++) { drtable[i] = 1.0/(rtable[i+1] - rtable[i]); dftable[i] = ftable[i+1] - ftable[i]; dctable[i] = ctable[i+1] - ctable[i]; detable[i] = etable[i+1] - etable[i]; } if (cut_respa) { for (int i = 0; i < ntablem1; i++) { dvtable[i] = vtable[i+1] - vtable[i]; dptable[i] = ptable[i+1] - ptable[i]; } } // get the delta values for the last table entries // tables are connected periodically between 0 and ntablem1 drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); dftable[ntablem1] = ftable[0] - ftable[ntablem1]; dctable[ntablem1] = ctable[0] - ctable[ntablem1]; detable[ntablem1] = etable[0] - etable[ntablem1]; if (cut_respa) { dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; dptable[ntablem1] = ptable[0] - ptable[ntablem1]; } // get the correct delta values at itablemax // smallest r is in bin itablemin // largest r is in bin itablemax, which is itablemin-1, // or ntablem1 if itablemin=0 // deltas at itablemax only needed if corresponding rsq < cut*cut // if so, compute deltas between rsq and cut*cut double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp; p_tmp = 0.0; v_tmp = 0.0; itablemin = minrsq_lookup.i & ncoulmask; itablemin >>= ncoulshiftbits; int itablemax = itablemin - 1; if (itablemin == 0) itablemax = ntablem1; rsq_lookup.i = itablemax << ncoulshiftbits; rsq_lookup.i |= maskhi; if (rsq_lookup.f < cut_coulsq) { rsq_lookup.f = cut_coulsq; r = sqrtf(rsq_lookup.f); if (msmflag) { egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)* force->kspace->dgamma(r/cut_coul); } else { grij = g_ewald * r; expm2 = exp(-grij*grij); derfc = erfc(grij); } if (cut_respa == NULL) { c_tmp = qqrd2e/r; if (msmflag) { f_tmp = qqrd2e/r * fgamma; e_tmp = qqrd2e/r * egamma; } else { f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); e_tmp = qqrd2e/r * derfc; } } else { c_tmp = 0.0; p_tmp = qqrd2e/r; if (msmflag) { f_tmp = qqrd2e/r * (fgamma - 1.0); e_tmp = qqrd2e/r * egamma; v_tmp = qqrd2e/r * fgamma; } else { f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2 - 1.0); e_tmp = qqrd2e/r * derfc; v_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); } if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) { if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) { rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]); f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw); } else { if (msmflag) f_tmp = qqrd2e/r * fgamma; else f_tmp = qqrd2e/r * (derfc + MY_ISPI4*grij*expm2); c_tmp = qqrd2e/r; } } } drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]); dftable[itablemax] = f_tmp - ftable[itablemax]; dctable[itablemax] = c_tmp - ctable[itablemax]; detable[itablemax] = e_tmp - etable[itablemax]; if (cut_respa) { dvtable[itablemax] = v_tmp - vtable[itablemax]; dptable[itablemax] = p_tmp - ptable[itablemax]; } } } /* ---------------------------------------------------------------------- setup force tables for dispersion used in compute routines ------------------------------------------------------------------------- */ void Pair::init_tables_disp(double cut_lj_global) { int masklo,maskhi; double rsq; double g_ewald_6 = force->kspace->g_ewald_6; double g2 = g_ewald_6*g_ewald_6, g6 = g2*g2*g2, g8 = g6*g2; tabinnerdispsq = tabinner_disp*tabinner_disp; init_bitmap(tabinner_disp,cut_lj_global,ndisptablebits, masklo,maskhi,ndispmask,ndispshiftbits); int ntable = 1; for (int i = 0; i < ndisptablebits; i++) ntable *= 2; // linear lookup tables of length N = 2^ndisptablebits // stored value = value at lower edge of bin // d values = delta from lower edge to upper edge of bin if (fdisptable) free_disp_tables(); memory->create(rdisptable,ntable,"pair:rdisptable"); memory->create(fdisptable,ntable,"pair:fdisptable"); memory->create(edisptable,ntable,"pair:edisptable"); memory->create(drdisptable,ntable,"pair:drdisptable"); memory->create(dfdisptable,ntable,"pair:dfdisptable"); memory->create(dedisptable,ntable,"pair:dedisptable"); union_int_float_t rsq_lookup; union_int_float_t minrsq_lookup; int itablemin; minrsq_lookup.i = 0 << ndispshiftbits; minrsq_lookup.i |= maskhi; for (int i = 0; i < ntable; i++) { rsq_lookup.i = i << ndispshiftbits; rsq_lookup.i |= masklo; if (rsq_lookup.f < tabinnerdispsq) { rsq_lookup.i = i << ndispshiftbits; rsq_lookup.i |= maskhi; } rsq = rsq_lookup.f; register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2); rdisptable[i] = rsq_lookup.f; fdisptable[i] = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; edisptable[i] = g6*((a2+1.0)*a2+0.5)*x2; minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); } tabinnerdispsq = minrsq_lookup.f; int ntablem1 = ntable - 1; for (int i = 0; i < ntablem1; i++) { drdisptable[i] = 1.0/(rdisptable[i+1] - rdisptable[i]); dfdisptable[i] = fdisptable[i+1] - fdisptable[i]; dedisptable[i] = edisptable[i+1] - edisptable[i]; } // get the delta values for the last table entries // tables are connected periodically between 0 and ntablem1 drdisptable[ntablem1] = 1.0/(rdisptable[0] - rdisptable[ntablem1]); dfdisptable[ntablem1] = fdisptable[0] - fdisptable[ntablem1]; dedisptable[ntablem1] = edisptable[0] - edisptable[ntablem1]; // get the correct delta values at itablemax // smallest r is in bin itablemin // largest r is in bin itablemax, which is itablemin-1, // or ntablem1 if itablemin=0 // deltas at itablemax only needed if corresponding rsq < cut*cut // if so, compute deltas between rsq and cut*cut double f_tmp,e_tmp; double cut_lj_globalsq; itablemin = minrsq_lookup.i & ndispmask; itablemin >>= ndispshiftbits; int itablemax = itablemin - 1; if (itablemin == 0) itablemax = ntablem1; rsq_lookup.i = itablemax << ndispshiftbits; rsq_lookup.i |= maskhi; if (rsq_lookup.f < (cut_lj_globalsq = cut_lj_global * cut_lj_global)) { rsq_lookup.f = cut_lj_globalsq; register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2); f_tmp = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; e_tmp = g6*((a2+1.0)*a2+0.5)*x2; drdisptable[itablemax] = 1.0/(rsq_lookup.f - rdisptable[itablemax]); dfdisptable[itablemax] = f_tmp - fdisptable[itablemax]; dedisptable[itablemax] = e_tmp - edisptable[itablemax]; } } /* ---------------------------------------------------------------------- free memory for tables used in Coulombic pair computations ------------------------------------------------------------------------- */ void Pair::free_tables() { memory->destroy(rtable); memory->destroy(drtable); memory->destroy(ftable); memory->destroy(dftable); memory->destroy(ctable); memory->destroy(dctable); memory->destroy(etable); memory->destroy(detable); memory->destroy(vtable); memory->destroy(dvtable); memory->destroy(ptable); memory->destroy(dptable); } /* ---------------------------------------------------------------------- free memory for tables used in pair computations for dispersion ------------------------------------------------------------------------- */ void Pair::free_disp_tables() { memory->destroy(rdisptable); memory->destroy(drdisptable); memory->destroy(fdisptable); memory->destroy(dfdisptable); memory->destroy(edisptable); memory->destroy(dedisptable); } /* ---------------------------------------------------------------------- mixing of pair potential prefactors (epsilon) ------------------------------------------------------------------------- */ double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2) { if (mix_flag == GEOMETRIC) return sqrt(eps1*eps2); else if (mix_flag == ARITHMETIC) return sqrt(eps1*eps2); else if (mix_flag == SIXTHPOWER) return (2.0 * sqrt(eps1*eps2) * pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0))); else return 0.0; } /* ---------------------------------------------------------------------- mixing of pair potential distances (sigma, cutoff) ------------------------------------------------------------------------- */ double Pair::mix_distance(double sig1, double sig2) { if (mix_flag == GEOMETRIC) return sqrt(sig1*sig2); else if (mix_flag == ARITHMETIC) return (0.5 * (sig1+sig2)); else if (mix_flag == SIXTHPOWER) return pow((0.5 * (pow(sig1,6.0) + pow(sig2,6.0))),1.0/6.0); else return 0.0; } /* ---------------------------------------------------------------------- */ void Pair::compute_dummy(int eflag, int vflag) { if (eflag || vflag) ev_setup(eflag,vflag); else evflag = 0; } /* ------------------------------------------------------------------- register a callback to a compute, so it can compute and accumulate additional properties during the pair computation from within Pair::ev_tally(). ensure each compute instance is registered only once ---------------------------------------------------------------------- */ void Pair::add_tally_callback(Compute *ptr) { if (lmp->kokkos) error->all(FLERR,"Cannot yet use compute tally with Kokkos"); int i,found=-1; for (i=0; i < num_tally_compute; ++i) { if (list_tally_compute[i] == ptr) found = i; } if (found < 0) { found = num_tally_compute; ++num_tally_compute; void *p = memory->srealloc((void *)list_tally_compute, sizeof(Compute *) * num_tally_compute, "pair:list_tally_compute"); list_tally_compute = (Compute **) p; list_tally_compute[num_tally_compute-1] = ptr; } } /* ------------------------------------------------------------------- unregister a callback to a fix for additional pairwise tallying ---------------------------------------------------------------------- */ void Pair::del_tally_callback(Compute *ptr) { int i,found=-1; for (i=0; i < num_tally_compute; ++i) { if (list_tally_compute[i] == ptr) found = i; } if (found < 0) return; // compact the list of active computes --num_tally_compute; for (i=found; i < num_tally_compute; ++i) { list_tally_compute[i] = list_tally_compute[i+1]; } } /* ---------------------------------------------------------------------- setup for energy, virial computation see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) ------------------------------------------------------------------------- */ void Pair::ev_setup(int eflag, int vflag) { int i,n; evflag = 1; eflag_either = eflag; eflag_global = eflag % 2; eflag_atom = eflag / 2; vflag_either = vflag; vflag_global = vflag % 4; vflag_atom = vflag / 4; // reallocate per-atom arrays if necessary if (eflag_atom && atom->nmax > maxeatom) { maxeatom = atom->nmax; memory->destroy(eatom); memory->create(eatom,comm->nthreads*maxeatom,"pair:eatom"); } if (vflag_atom && atom->nmax > maxvatom) { maxvatom = atom->nmax; memory->destroy(vatom); memory->create(vatom,comm->nthreads*maxvatom,6,"pair:vatom"); } // zero accumulators // use force->newton instead of newton_pair // b/c some bonds/dihedrals call pair::ev_tally with pairwise info if (eflag_global) eng_vdwl = eng_coul = 0.0; if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; if (eflag_atom) { n = atom->nlocal; if (force->newton) n += atom->nghost; for (i = 0; i < n; i++) eatom[i] = 0.0; } if (vflag_atom) { n = atom->nlocal; if (force->newton) n += atom->nghost; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; vatom[i][1] = 0.0; vatom[i][2] = 0.0; vatom[i][3] = 0.0; vatom[i][4] = 0.0; vatom[i][5] = 0.0; } } // if vflag_global = 2 and pair::compute() calls virial_fdotr_compute() // compute global virial via (F dot r) instead of via pairwise summation // unset other flags as appropriate if (vflag_global == 2 && no_virial_fdotr_compute == 0) { vflag_fdotr = 1; vflag_global = 0; if (vflag_atom == 0) vflag_either = 0; if (vflag_either == 0 && eflag_either == 0) evflag = 0; } else vflag_fdotr = 0; } /* ---------------------------------------------------------------------- set all flags to zero for energy, virial computation called by some complicated many-body potentials that use individual flags to insure no holdover of flags from previous timestep ------------------------------------------------------------------------- */ void Pair::ev_unset() { evflag = 0; eflag_either = 0; eflag_global = 0; eflag_atom = 0; vflag_either = 0; vflag_global = 0; vflag_atom = 0; vflag_fdotr = 0; } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators need i < nlocal test since called by bond_quartic and dihedral_charmm ------------------------------------------------------------------------- */ void Pair::ev_tally(int i, int j, int nlocal, int newton_pair, double evdwl, double ecoul, double fpair, double delx, double dely, double delz) { double evdwlhalf,ecoulhalf,epairhalf,v[6]; if (eflag_either) { if (eflag_global) { if (newton_pair) { eng_vdwl += evdwl; eng_coul += ecoul; } else { evdwlhalf = 0.5*evdwl; ecoulhalf = 0.5*ecoul; if (i < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } if (j < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } } } if (eflag_atom) { epairhalf = 0.5 * (evdwl + ecoul); if (newton_pair || i < nlocal) eatom[i] += epairhalf; if (newton_pair || j < nlocal) eatom[j] += epairhalf; } } if (vflag_either) { v[0] = delx*delx*fpair; v[1] = dely*dely*fpair; v[2] = delz*delz*fpair; v[3] = delx*dely*fpair; v[4] = delx*delz*fpair; v[5] = dely*delz*fpair; if (vflag_global) { if (newton_pair) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } else { if (i < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } if (j < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } } } if (vflag_atom) { if (newton_pair || i < nlocal) { vatom[i][0] += 0.5*v[0]; vatom[i][1] += 0.5*v[1]; vatom[i][2] += 0.5*v[2]; vatom[i][3] += 0.5*v[3]; vatom[i][4] += 0.5*v[4]; vatom[i][5] += 0.5*v[5]; } if (newton_pair || j < nlocal) { vatom[j][0] += 0.5*v[0]; vatom[j][1] += 0.5*v[1]; vatom[j][2] += 0.5*v[2]; vatom[j][3] += 0.5*v[3]; vatom[j][4] += 0.5*v[4]; vatom[j][5] += 0.5*v[5]; } } } if (num_tally_compute > 0) { for (int k=0; k < num_tally_compute; ++k) { Compute *c = list_tally_compute[k]; c->pair_tally_callback(i, j, nlocal, newton_pair, evdwl, ecoul, fpair, delx, dely, delz); } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators can use this version with full neighbor lists ------------------------------------------------------------------------- */ void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair, double delx, double dely, double delz) { double v[6]; if (eflag_either) { if (eflag_global) { eng_vdwl += 0.5*evdwl; eng_coul += 0.5*ecoul; } if (eflag_atom) eatom[i] += 0.5 * (evdwl + ecoul); } if (vflag_either) { v[0] = 0.5*delx*delx*fpair; v[1] = 0.5*dely*dely*fpair; v[2] = 0.5*delz*delz*fpair; v[3] = 0.5*delx*dely*fpair; v[4] = 0.5*delx*delz*fpair; v[5] = 0.5*dely*delz*fpair; if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz ------------------------------------------------------------------------- */ void Pair::ev_tally_xyz(int i, int j, int nlocal, int newton_pair, double evdwl, double ecoul, double fx, double fy, double fz, double delx, double dely, double delz) { double evdwlhalf,ecoulhalf,epairhalf,v[6]; if (eflag_either) { if (eflag_global) { if (newton_pair) { eng_vdwl += evdwl; eng_coul += ecoul; } else { evdwlhalf = 0.5*evdwl; ecoulhalf = 0.5*ecoul; if (i < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } if (j < nlocal) { eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } } } if (eflag_atom) { epairhalf = 0.5 * (evdwl + ecoul); if (newton_pair || i < nlocal) eatom[i] += epairhalf; if (newton_pair || j < nlocal) eatom[j] += epairhalf; } } if (vflag_either) { v[0] = delx*fx; v[1] = dely*fy; v[2] = delz*fz; v[3] = delx*fy; v[4] = delx*fz; v[5] = dely*fz; if (vflag_global) { if (newton_pair) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } else { if (i < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } if (j < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } } } if (vflag_atom) { if (newton_pair || i < nlocal) { vatom[i][0] += 0.5*v[0]; vatom[i][1] += 0.5*v[1]; vatom[i][2] += 0.5*v[2]; vatom[i][3] += 0.5*v[3]; vatom[i][4] += 0.5*v[4]; vatom[i][5] += 0.5*v[5]; } if (newton_pair || j < nlocal) { vatom[j][0] += 0.5*v[0]; vatom[j][1] += 0.5*v[1]; vatom[j][2] += 0.5*v[2]; vatom[j][3] += 0.5*v[3]; vatom[j][4] += 0.5*v[4]; vatom[j][5] += 0.5*v[5]; } } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators for virial, have delx,dely,delz and fx,fy,fz called when using full neighbor lists ------------------------------------------------------------------------- */ void Pair::ev_tally_xyz_full(int i, double evdwl, double ecoul, double fx, double fy, double fz, double delx, double dely, double delz) { double evdwlhalf,ecoulhalf,epairhalf,v[6]; if (eflag_either) { if (eflag_global) { evdwlhalf = 0.5*evdwl; ecoulhalf = 0.5*ecoul; eng_vdwl += evdwlhalf; eng_coul += ecoulhalf; } if (eflag_atom) { epairhalf = 0.5 * (evdwl + ecoul); eatom[i] += epairhalf; } } if (vflag_either) { v[0] = 0.5*delx*fx; v[1] = 0.5*dely*fy; v[2] = 0.5*delz*fz; v[3] = 0.5*delx*fy; v[4] = 0.5*delx*fz; v[5] = 0.5*dely*fz; if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators called by SW and hbond potentials, newton_pair is always on virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk ------------------------------------------------------------------------- */ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, double *fj, double *fk, double *drji, double *drki) { double epairthird,v[6]; if (eflag_either) { if (eflag_global) { eng_vdwl += evdwl; eng_coul += ecoul; } if (eflag_atom) { epairthird = THIRD * (evdwl + ecoul); eatom[i] += epairthird; eatom[j] += epairthird; eatom[k] += epairthird; } } if (vflag_either) { v[0] = drji[0]*fj[0] + drki[0]*fk[0]; v[1] = drji[1]*fj[1] + drki[1]*fk[1]; v[2] = drji[2]*fj[2] + drki[2]*fk[2]; v[3] = drji[0]*fj[1] + drki[0]*fk[1]; v[4] = drji[0]*fj[2] + drki[0]*fk[2]; v[5] = drji[1]*fj[2] + drki[1]*fk[2]; if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { vatom[i][0] += THIRD*v[0]; vatom[i][1] += THIRD*v[1]; vatom[i][2] += THIRD*v[2]; vatom[i][3] += THIRD*v[3]; vatom[i][4] += THIRD*v[4]; vatom[i][5] += THIRD*v[5]; vatom[j][0] += THIRD*v[0]; vatom[j][1] += THIRD*v[1]; vatom[j][2] += THIRD*v[2]; vatom[j][3] += THIRD*v[3]; vatom[j][4] += THIRD*v[4]; vatom[j][5] += THIRD*v[5]; vatom[k][0] += THIRD*v[0]; vatom[k][1] += THIRD*v[1]; vatom[k][2] += THIRD*v[2]; vatom[k][3] += THIRD*v[3]; vatom[k][4] += THIRD*v[4]; vatom[k][5] += THIRD*v[5]; } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::ev_tally4(int i, int j, int k, int m, double evdwl, double *fi, double *fj, double *fk, double *drim, double *drjm, double *drkm) { double epairfourth,v[6]; if (eflag_either) { if (eflag_global) eng_vdwl += evdwl; if (eflag_atom) { epairfourth = 0.25 * evdwl; eatom[i] += epairfourth; eatom[j] += epairfourth; eatom[k] += epairfourth; eatom[m] += epairfourth; } } if (vflag_atom) { v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; } } /* ---------------------------------------------------------------------- tally ecoul and virial into each of atoms in list called by TIP4P potential, newton_pair is always on weight assignments by alpha, so contribution is all to O atom as alpha -> 0.0 key = 0 if neither atom = water O key = 1 if first atom = water O key = 2 if second atom = water O key = 3 if both atoms = water O ------------------------------------------------------------------------- */ void Pair::ev_tally_tip4p(int key, int *list, double *v, double ecoul, double alpha) { int i; if (eflag_either) { if (eflag_global) eng_coul += ecoul; if (eflag_atom) { if (key == 0) { eatom[list[0]] += 0.5*ecoul; eatom[list[1]] += 0.5*ecoul; } else if (key == 1) { eatom[list[0]] += 0.5*ecoul*(1-alpha); eatom[list[1]] += 0.25*ecoul*alpha; eatom[list[2]] += 0.25*ecoul*alpha; eatom[list[3]] += 0.5*ecoul; } else if (key == 2) { eatom[list[0]] += 0.5*ecoul; eatom[list[1]] += 0.5*ecoul*(1-alpha); eatom[list[2]] += 0.25*ecoul*alpha; eatom[list[3]] += 0.25*ecoul*alpha; } else { eatom[list[0]] += 0.5*ecoul*(1-alpha); eatom[list[1]] += 0.25*ecoul*alpha; eatom[list[2]] += 0.25*ecoul*alpha; eatom[list[3]] += 0.5*ecoul*(1-alpha); eatom[list[4]] += 0.25*ecoul*alpha; eatom[list[5]] += 0.25*ecoul*alpha; } } } if (vflag_either) { if (vflag_global) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } if (vflag_atom) { if (key == 0) { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]; vatom[list[1]][i] += 0.5*v[i]; } } else if (key == 1) { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]*(1-alpha); vatom[list[1]][i] += 0.25*v[i]*alpha; vatom[list[2]][i] += 0.25*v[i]*alpha; vatom[list[3]][i] += 0.5*v[i]; } } else if (key == 2) { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]; vatom[list[1]][i] += 0.5*v[i]*(1-alpha); vatom[list[2]][i] += 0.25*v[i]*alpha; vatom[list[3]][i] += 0.25*v[i]*alpha; } } else { for (i = 0; i <= 5; i++) { vatom[list[0]][i] += 0.5*v[i]*(1-alpha); vatom[list[1]][i] += 0.25*v[i]*alpha; vatom[list[2]][i] += 0.25*v[i]*alpha; vatom[list[3]][i] += 0.5*v[i]*(1-alpha); vatom[list[4]][i] += 0.25*v[i]*alpha; vatom[list[5]][i] += 0.25*v[i]*alpha; } } } } } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by REAX/C potential, newton_pair is always on fi is magnitude of force on atom i ------------------------------------------------------------------------- */ void Pair::v_tally(int i, double *fi, double *deli) { double v[6]; v[0] = 0.5*deli[0]*fi[0]; v[1] = 0.5*deli[1]*fi[1]; v[2] = 0.5*deli[2]*fi[2]; v[3] = 0.5*deli[0]*fi[1]; v[4] = 0.5*deli[0]*fi[2]; v[5] = 0.5*deli[1]*fi[2]; vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by AIREBO potential, newton_pair is always on fpair is magnitude of force on atom I ------------------------------------------------------------------------- */ void Pair::v_tally2(int i, int j, double fpair, double *drij) { double v[6]; v[0] = 0.5 * drij[0]*drij[0]*fpair; v[1] = 0.5 * drij[1]*drij[1]*fpair; v[2] = 0.5 * drij[2]*drij[2]*fpair; v[3] = 0.5 * drij[0]*drij[1]*fpair; v[4] = 0.5 * drij[0]*drij[2]*fpair; v[5] = 0.5 * drij[1]*drij[2]*fpair; vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by AIREBO and Tersoff potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally3(int i, int j, int k, double *fi, double *fj, double *drik, double *drjk) { double v[6]; v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into per-atom accumulators called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally4(int i, int j, int k, int m, double *fi, double *fj, double *fk, double *drim, double *drjm, double *drkm) { double v[6]; v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; } /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators called by pair lubricate potential with 6 tensor components ------------------------------------------------------------------------- */ void Pair::v_tally_tensor(int i, int j, int nlocal, int newton_pair, double vxx, double vyy, double vzz, double vxy, double vxz, double vyz) { double v[6]; v[0] = vxx; v[1] = vyy; v[2] = vzz; v[3] = vxy; v[4] = vxz; v[5] = vyz; if (vflag_global) { if (newton_pair) { virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2]; virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5]; } else { if (i < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } if (j < nlocal) { virial[0] += 0.5*v[0]; virial[1] += 0.5*v[1]; virial[2] += 0.5*v[2]; virial[3] += 0.5*v[3]; virial[4] += 0.5*v[4]; virial[5] += 0.5*v[5]; } } } if (vflag_atom) { if (newton_pair || i < nlocal) { vatom[i][0] += 0.5*v[0]; vatom[i][1] += 0.5*v[1]; vatom[i][2] += 0.5*v[2]; vatom[i][3] += 0.5*v[3]; vatom[i][4] += 0.5*v[4]; vatom[i][5] += 0.5*v[5]; } if (newton_pair || j < nlocal) { vatom[j][0] += 0.5*v[0]; vatom[j][1] += 0.5*v[1]; vatom[j][2] += 0.5*v[2]; vatom[j][3] += 0.5*v[3]; vatom[j][4] += 0.5*v[4]; vatom[j][5] += 0.5*v[5]; } } } /* ---------------------------------------------------------------------- compute global pair virial via summing F dot r over own & ghost atoms at this point, only pairwise forces have been accumulated in atom->f ------------------------------------------------------------------------- */ void Pair::virial_fdotr_compute() { double **x = atom->x; double **f = atom->f; // sum over force on all particles including ghosts if (neighbor->includegroup == 0) { int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) { virial[0] += f[i][0]*x[i][0]; virial[1] += f[i][1]*x[i][1]; virial[2] += f[i][2]*x[i][2]; virial[3] += f[i][1]*x[i][0]; virial[4] += f[i][2]*x[i][0]; virial[5] += f[i][2]*x[i][1]; } // neighbor includegroup flag is set // sum over force on initial nfirst particles and ghosts } else { int nall = atom->nfirst; for (int i = 0; i < nall; i++) { virial[0] += f[i][0]*x[i][0]; virial[1] += f[i][1]*x[i][1]; virial[2] += f[i][2]*x[i][2]; virial[3] += f[i][1]*x[i][0]; virial[4] += f[i][2]*x[i][0]; virial[5] += f[i][2]*x[i][1]; } nall = atom->nlocal + atom->nghost; for (int i = atom->nlocal; i < nall; i++) { virial[0] += f[i][0]*x[i][0]; virial[1] += f[i][1]*x[i][1]; virial[2] += f[i][2]*x[i][2]; virial[3] += f[i][1]*x[i][0]; virial[4] += f[i][2]*x[i][0]; virial[5] += f[i][2]*x[i][1]; } } // prevent multiple calls to update the virial // when a hybrid pair style uses both a gpu and non-gpu pair style // or when respa is used with gpu pair styles vflag_fdotr = 0; } /* ---------------------------------------------------------------------- write a table of pair potential energy/force vs distance to a file ------------------------------------------------------------------------- */ void Pair::write_file(int narg, char **arg) { if (narg < 8) error->all(FLERR,"Illegal pair_write command"); if (single_enable == 0) error->all(FLERR,"Pair style does not support pair_write"); // parse arguments int itype = force->inumeric(FLERR,arg[0]); int jtype = force->inumeric(FLERR,arg[1]); if (itype < 1 || itype > atom->ntypes || jtype < 1 || jtype > atom->ntypes) error->all(FLERR,"Invalid atom types in pair_write command"); int n = force->inumeric(FLERR,arg[2]); int style = NONE; if (strcmp(arg[3],"r") == 0) style = RLINEAR; else if (strcmp(arg[3],"rsq") == 0) style = RSQ; else if (strcmp(arg[3],"bitmap") == 0) style = BMP; else error->all(FLERR,"Invalid style in pair_write command"); double inner = force->numeric(FLERR,arg[4]); double outer = force->numeric(FLERR,arg[5]); if (inner <= 0.0 || inner >= outer) error->all(FLERR,"Invalid cutoffs in pair_write command"); // open file in append mode // print header in format used by pair_style table int me; MPI_Comm_rank(world,&me); FILE *fp; if (me == 0) { fp = fopen(arg[6],"a"); if (fp == NULL) error->one(FLERR,"Cannot open pair_write file"); fprintf(fp,"# Pair potential %s for atom types %d %d: i,r,energy,force\n", force->pair_style,itype,jtype); if (style == RLINEAR) fprintf(fp,"\n%s\nN %d R %.15g %.15g\n\n",arg[7],n,inner,outer); if (style == RSQ) fprintf(fp,"\n%s\nN %d RSQ %.15g %.15g\n\n",arg[7],n,inner,outer); } // initialize potentials before evaluating pair potential // insures all pair coeffs are set and force constants // also initialize neighbor so that neighbor requests are processed // NOTE: might be safest to just do lmp->init() force->init(); neighbor->init(); // if pair style = any of EAM, swap in dummy fp vector double eamfp[2]; eamfp[0] = eamfp[1] = 0.0; double *eamfp_hold; Pair *epair = force->pair_match("eam",0); if (epair) epair->swap_eam(eamfp,&eamfp_hold); // if atom style defines charge, swap in dummy q vec double q[2]; q[0] = q[1] = 1.0; if (narg == 10) { q[0] = force->numeric(FLERR,arg[8]); q[1] = force->numeric(FLERR,arg[9]); } double *q_hold; if (atom->q) { q_hold = atom->q; atom->q = q; } // evaluate energy and force at each of N distances int masklo,maskhi,nmask,nshiftbits; if (style == BMP) { init_bitmap(inner,outer,n,masklo,maskhi,nmask,nshiftbits); int ntable = 1 << n; if (me == 0) fprintf(fp,"\n%s\nN %d BITMAP %.15g %.15g\n\n",arg[7],ntable,inner,outer); n = ntable; } double r,e,f,rsq; union_int_float_t rsq_lookup; for (int i = 0; i < n; i++) { if (style == RLINEAR) { r = inner + (outer-inner) * i/(n-1); rsq = r*r; } else if (style == RSQ) { rsq = inner*inner + (outer*outer - inner*inner) * i/(n-1); r = sqrt(rsq); } else if (style == BMP) { rsq_lookup.i = i << nshiftbits; rsq_lookup.i |= masklo; if (rsq_lookup.f < inner*inner) { rsq_lookup.i = i << nshiftbits; rsq_lookup.i |= maskhi; } rsq = rsq_lookup.f; r = sqrt(rsq); } if (rsq < cutsq[itype][jtype]) { e = single(0,1,itype,jtype,rsq,1.0,1.0,f); f *= r; } else e = f = 0.0; if (me == 0) fprintf(fp,"%d %.15g %.15g %.15g\n",i+1,r,e,f); } // restore original vecs that were swapped in for double *tmp; if (epair) epair->swap_eam(eamfp_hold,&tmp); if (atom->q) atom->q = q_hold; if (me == 0) fclose(fp); } /* ---------------------------------------------------------------------- define bitmap parameters based on inner and outer cutoffs ------------------------------------------------------------------------- */ void Pair::init_bitmap(double inner, double outer, int ntablebits, int &masklo, int &maskhi, int &nmask, int &nshiftbits) { if (sizeof(int) != sizeof(float)) error->all(FLERR,"Bitmapped lookup tables require int/float be same size"); if (ntablebits > sizeof(float)*CHAR_BIT) error->all(FLERR,"Too many total bits for bitmapped lookup table"); if (inner >= outer) error->warning(FLERR,"Table inner cutoff >= outer cutoff"); int nlowermin = 1; while (!((pow(double(2),(double)nlowermin) <= inner*inner) && (pow(double(2),(double)nlowermin+1.0) > inner*inner))) { if (pow(double(2),(double)nlowermin) <= inner*inner) nlowermin++; else nlowermin--; } int nexpbits = 0; double required_range = outer*outer / pow(double(2),(double)nlowermin); double available_range = 2.0; while (available_range < required_range) { nexpbits++; available_range = pow(double(2),pow(double(2),(double)nexpbits)); } int nmantbits = ntablebits - nexpbits; if (nexpbits > sizeof(float)*CHAR_BIT - FLT_MANT_DIG) error->all(FLERR,"Too many exponent bits for lookup table"); if (nmantbits+1 > FLT_MANT_DIG) error->all(FLERR,"Too many mantissa bits for lookup table"); if (nmantbits < 3) error->all(FLERR,"Too few bits for lookup table"); nshiftbits = FLT_MANT_DIG - (nmantbits+1); nmask = 1; for (int j = 0; j < ntablebits+nshiftbits; j++) nmask *= 2; nmask -= 1; union_int_float_t rsq_lookup; rsq_lookup.f = outer*outer; maskhi = rsq_lookup.i & ~(nmask); rsq_lookup.f = inner*inner; masklo = rsq_lookup.i & ~(nmask); } /* ---------------------------------------------------------------------- */ double Pair::memory_usage() { double bytes = comm->nthreads*maxeatom * sizeof(double); bytes += comm->nthreads*maxvatom*6 * sizeof(double); return bytes; }