From ba728e5e07e11914af881eef2ae829e6d41f9fcd Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 1 Dec 2011 15:47:17 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7262 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/KSPACE/ewald.cpp | 20 ++++++++------- src/KSPACE/ewald.h | 1 - src/KSPACE/pppm.cpp | 16 ++++++------ src/KSPACE/pppm.h | 1 - src/KSPACE/pppm_cg.cpp | 16 +++++++----- src/KSPACE/pppm_tip4p.cpp | 3 ++- src/MANYBODY/pair_airebo.cpp | 4 +-- src/MANYBODY/pair_airebo.h | 2 +- src/MANYBODY/pair_comb.cpp | 33 ++++-------------------- src/MANYBODY/pair_comb.h | 46 ++++++++++++++++++++++++++-------- src/MOLECULE/angle_hybrid.cpp | 10 ++++++++ src/MOLECULE/angle_hybrid.h | 8 +++--- src/MOLECULE/dihedral_hybrid.h | 7 +++--- src/MOLECULE/improper_hybrid.h | 7 +++--- src/USER-EWALDN/ewald_n.cpp | 22 ++++++++++------ src/USER-EWALDN/ewald_n.h | 2 +- src/pair.h | 2 ++ 17 files changed, 115 insertions(+), 85 deletions(-) diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index b6cff3d932..4e2782ff04 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -95,7 +95,6 @@ void Ewald::init() // extract short-range Coulombic cutoff from pair style - qqrd2e = force->qqrd2e; scale = 1.0; if (force->pair == NULL) @@ -270,10 +269,12 @@ void Ewald::compute(int eflag, int vflag) // convert E-field to force + const double qscale = force->qqrd2e * scale; + for (i = 0; i < nlocal; i++) { - f[i][0] += qqrd2e*scale * q[i]*ek[i][0]; - f[i][1] += qqrd2e*scale * q[i]*ek[i][1]; - f[i][2] += qqrd2e*scale * q[i]*ek[i][2]; + f[i][0] += qscale * q[i]*ek[i][0]; + f[i][1] += qscale * q[i]*ek[i][1]; + f[i][2] += qscale * q[i]*ek[i][2]; } // energy if requested @@ -284,7 +285,7 @@ void Ewald::compute(int eflag, int vflag) sfacim_all[k]*sfacim_all[k]); energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); - energy *= qqrd2e*scale; + energy *= qscale; } // virial if requested @@ -295,7 +296,7 @@ void Ewald::compute(int eflag, int vflag) uk = ug[k] * (sfacrl_all[k]*sfacrl_all[k] + sfacim_all[k]*sfacim_all[k]); for (n = 0; n < 6; n++) virial[n] += uk*vg[k][n]; } - for (n = 0; n < 6; n++) virial[n] *= qqrd2e*scale; + for (n = 0; n < 6; n++) virial[n] *= qscale; } if (slabflag) slabcorr(eflag); @@ -817,16 +818,17 @@ void Ewald::slabcorr(int eflag) // compute corrections - double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double qscale = force->qqrd2e * scale; - if (eflag) energy += qqrd2e*scale * e_slabcorr; + if (eflag) energy += qscale * e_slabcorr; // add on force corrections double ffact = -4.0*MY_PI*dipole_all/volume; double **f = atom->f; - for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact; + for (int i = 0; i < nlocal; i++) f[i][2] += qscale * q[i]*ffact; } /* ---------------------------------------------------------------------- diff --git a/src/KSPACE/ewald.h b/src/KSPACE/ewald.h index cd22f987f7..e88b8989ac 100644 --- a/src/KSPACE/ewald.h +++ b/src/KSPACE/ewald.h @@ -36,7 +36,6 @@ class Ewald : public KSpace { protected: double precision; int kcount,kmax,kmax3d,kmax_created; - double qqrd2e; double gsqmx,qsum,qsqsum,volume; int nmax; diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index e2d017a6e5..4e810bbc23 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -139,7 +139,6 @@ void PPPM::init() // extract short-range Coulombic cutoff from pair style - qqrd2e = force->qqrd2e; scale = 1.0; if (force->pair == NULL) @@ -714,6 +713,8 @@ void PPPM::compute(int eflag, int vflag) // sum energy across procs and add in volume-dependent term + const double qscale = force->qqrd2e * scale; + if (eflag) { double energy_all; MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); @@ -722,7 +723,7 @@ void PPPM::compute(int eflag, int vflag) energy *= 0.5*volume; energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); - energy *= qqrd2e*scale; + energy *= qscale; } // sum virial across procs @@ -730,7 +731,7 @@ void PPPM::compute(int eflag, int vflag) if (vflag) { double virial_all[6]; MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i]; + for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i]; } // 2d slab correction @@ -1734,7 +1735,7 @@ void PPPM::fieldforce() } // convert E-field to force - const double qfactor = qqrd2e*scale*q[i]; + const double qfactor = force->qqrd2e * scale * q[i]; f[i][0] += qfactor*ekx; f[i][1] += qfactor*eky; f[i][2] += qfactor*ekz; @@ -1887,16 +1888,17 @@ void PPPM::slabcorr(int eflag) // compute corrections - double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double qscale = force->qqrd2e * scale; - if (eflag) energy += qqrd2e*scale * e_slabcorr; + if (eflag) energy += qscale * e_slabcorr; // add on force corrections double ffact = -4.0*MY_PI*dipole_all/volume; double **f = atom->f; - for (int i = 0; i < nlocal; i++) f[i][2] += qqrd2e*scale * q[i]*ffact; + for (int i = 0; i < nlocal; i++) f[i][2] += qscale * q[i]*ffact; } /* ---------------------------------------------------------------------- diff --git a/src/KSPACE/pppm.h b/src/KSPACE/pppm.h index e2a58e8804..1ee2cfbb4f 100644 --- a/src/KSPACE/pppm.h +++ b/src/KSPACE/pppm.h @@ -51,7 +51,6 @@ class PPPM : public KSpace { int nfactors; int *factors; double qsum,qsqsum; - double qqrd2e; double cutoff; double volume; double delxinv,delyinv,delzinv,delvolinv; diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp index daffb5ff2d..663d263a20 100644 --- a/src/KSPACE/pppm_cg.cpp +++ b/src/KSPACE/pppm_cg.cpp @@ -23,6 +23,7 @@ #include "atom.h" #include "domain.h" #include "error.h" +#include "force.h" #include "memory.h" #include "pppm_cg.h" @@ -172,6 +173,8 @@ void PPPMCG::compute(int eflag, int vflag) // sum energy across procs and add in volume-dependent term + const double qscale = force->qqrd2e * scale; + if (eflag) { double energy_all; MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world); @@ -180,7 +183,7 @@ void PPPMCG::compute(int eflag, int vflag) energy *= 0.5*volume; energy -= g_ewald*qsqsum/MY_PIS + MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume); - energy *= qqrd2e*scale; + energy *= qscale; } // sum virial across procs @@ -188,7 +191,7 @@ void PPPMCG::compute(int eflag, int vflag) if (vflag) { double virial_all[6]; MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world); - for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i]; + for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i]; } // 2d slab correction @@ -341,7 +344,7 @@ void PPPMCG::fieldforce() } // convert E-field to force - const double qfactor = qqrd2e*scale*q[i]; + const double qfactor = force->qqrd2e * scale * q[i]; f[i][0] += qfactor*ekx; f[i][1] += qfactor*eky; f[i][2] += qfactor*ekz; @@ -375,13 +378,14 @@ void PPPMCG::slabcorr(int eflag) // compute corrections - double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double e_slabcorr = 2.0*MY_PI*dipole_all*dipole_all/volume; + const double qscale = force->qqrd2e * scale; - if (eflag) energy += qqrd2e*scale * e_slabcorr; + if (eflag) energy += qscale * e_slabcorr; // add on force corrections - double ffact = -4.0*MY_PI*dipole_all/volume * qqrd2e * scale; + const double ffact = -4.0*MY_PI*dipole_all/volume * qscale; double **f = atom->f; for (int j = 0; j < num_charged; j++) { diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp index a191fa742c..16a0f7a04c 100644 --- a/src/KSPACE/pppm_tip4p.cpp +++ b/src/KSPACE/pppm_tip4p.cpp @@ -19,6 +19,7 @@ #include "pppm_tip4p.h" #include "atom.h" #include "domain.h" +#include "force.h" #include "memory.h" #include "error.h" @@ -205,7 +206,7 @@ void PPPMTIP4P::fieldforce() } // convert E-field to force - const double qfactor = qqrd2e*scale*q[i]; + const double qfactor = force->qqrd2e * scale * q[i]; if (type[i] != typeO) { f[i][0] += qfactor*ekx; f[i][1] += qfactor*eky; diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index c2625072cc..fcbbc1e091 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -3283,10 +3283,10 @@ double PairAIREBO::TijSpline(double Nij, double Nji, add pages to REBO neighbor list ------------------------------------------------------------------------- */ -void PairAIREBO::add_pages() +void PairAIREBO::add_pages(int howmany) { int toppage = maxpage; - maxpage += PGDELTA; + maxpage += howmany*PGDELTA; pages = (int **) memory->srealloc(pages,maxpage*sizeof(int *),"AIREBO:pages"); diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index 0c4795ab0d..5ef5d0ee14 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -101,7 +101,7 @@ class PairAIREBO : public Pair { double piRCSpline(double, double, double, int, int, double *); double TijSpline(double, double, double, double *); - void add_pages(); + void add_pages(int howmany=1); void read_file(char *); double Sp5th(double, double *, double *); diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index 51bc45ad93..68555bd4f4 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -1243,27 +1243,6 @@ double PairComb::comb_bij_d(double zeta, Param *param) /* ---------------------------------------------------------------------- */ -double PairComb::comb_gijk(double costheta, Param *param) -{ - double comb_c = param->c; - double comb_d = param->d; - - return (1.0 + pow(comb_c/comb_d,2.0) - - pow(comb_c,2.0) / (pow(comb_d,2.0) + pow(param->h - costheta,2.0))); -} - -/* ---------------------------------------------------------------------- */ - -double PairComb::comb_gijk_d(double costheta, Param *param) -{ - double numerator = -2.0 * pow(param->c,2) * (param->h - costheta); - double denominator = pow(pow(param->d,2.0) + - pow(param->h - costheta,2.0),2.0); - return numerator/denominator; -} - -/*------------------------------------------------------------------------- */ - void PairComb::attractive(Param *param, double prefactor, double rsqij, double rsqik, double *delrij, double *delrik, @@ -1649,10 +1628,10 @@ void PairComb::field(Param *param, double rsq, double iq,double jq, double PairComb::yasu_char(double *qf_fix, int &igroup) { - int i,j,k,ii,jj,kk,jnum,itag,jtag; - int itype,jtype,ktype,iparam_i,iparam_ij,iparam_ijk; + int i,j,ii,jj,jnum,itag,jtag; + int itype,jtype,iparam_i,iparam_ij; double xtmp,ytmp,ztmp; - double rsq1,rsq2,delr1[3],delr2[3],zeta_ij; + double rsq1,delr1[3]; int *ilist,*jlist,*numneigh,**firstneigh; double iq,jq,fqi,fqj,fqij,fqjj; double potal,fac11,fac11e,sr1,sr2,sr3; @@ -2147,15 +2126,13 @@ void PairComb::Short_neigh() /* ---------------------------------------------------------------------- */ -void PairComb::add_pages() +void PairComb::add_pages(int howmany) { int toppage = maxpage; - maxpage += PGDELTA; + maxpage += howmany*PGDELTA; pages = (int **) memory->srealloc(pages,maxpage*sizeof(int *),"pair:pages"); for (int i = toppage; i < maxpage; i++) memory->create(pages[i],pgsize,"pair:pages[i]"); } - -/* ---------------------------------------------------------------------- */ diff --git a/src/MANYBODY/pair_comb.h b/src/MANYBODY/pair_comb.h index 0a8b11ad05..19bc849b44 100644 --- a/src/MANYBODY/pair_comb.h +++ b/src/MANYBODY/pair_comb.h @@ -28,16 +28,16 @@ class PairComb : public Pair { public: PairComb(class LAMMPS *); virtual ~PairComb(); - void compute(int, int); + virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); void init_style(); double init_one(int, int); double memory_usage(); - double yasu_char(double *, int &); + virtual double yasu_char(double *, int &); - private: + protected: struct Param { double lam11,lam12,lam21,lam22; double c,d,h; @@ -100,8 +100,24 @@ class PairComb : public Pair { virtual double comb_fa_d(double, Param *, double,double); double comb_bij(double, Param *); double comb_bij_d(double, Param *); - double comb_gijk(double, Param *); - double comb_gijk_d(double, Param *); + + inline double comb_gijk(const double costheta, const Param * const param) const { + const double comb_c = param->c * param->c; + const double comb_d = param->d * param->d; + const double hcth = param->h - costheta; + + return param->gamma*(1.0 + comb_c/comb_d - comb_c / (comb_d + hcth*hcth)); + } + + inline double comb_gijk_d(const double costheta, const Param * const param) const { + const double comb_c = param->c * param->c; + const double comb_d = param->d * param->d; + const double hcth = param->h - costheta; + const double numerator = -2.0 * comb_c * hcth; + const double denominator = 1.0/(comb_d + hcth*hcth); + return param->gamma*numerator*denominator*denominator; + } + void comb_zetaterm_d(double, double *, double, double *, double, double *, double *, double *, Param *); void costheta_d(double *, double, double *, double, @@ -129,7 +145,7 @@ class PairComb : public Pair { // Short range neighbor list - void add_pages(); + void add_pages(int howmany=1); void Short_neigh(); int maxpage, pgsize, oneatom, **pages; int *sht_num, **sht_first; // short-range neighbor list @@ -137,17 +153,25 @@ class PairComb : public Pair { // vector functions, inline for efficiency - inline double vec3_dot(double *x, double *y) { + inline double vec3_dot(const double x[3], const double y[3]) const { return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; } - inline void vec3_add(double *x, double *y, double *z) { + + inline void vec3_add(const double x[3], const double y[3], + double * const z) const { z[0] = x[0]+y[0]; z[1] = x[1]+y[1]; z[2] = x[2]+y[2]; } - inline void vec3_scale(double k, double *x, double *y) { + + inline void vec3_scale(const double k, const double x[3], + double y[3]) const { y[0] = k*x[0]; y[1] = k*x[1]; y[2] = k*x[2]; } - inline void vec3_scaleadd(double k, double *x, double *y, double *z) { - z[0] = k*x[0]+y[0]; z[1] = k*x[1]+y[1]; z[2] = k*x[2]+y[2]; + + inline void vec3_scaleadd(const double k, const double x[3], + const double y[3], double * const z) const { + z[0] = k*x[0]+y[0]; + z[1] = k*x[1]+y[1]; + z[2] = k*x[2]+y[2]; } }; diff --git a/src/MOLECULE/angle_hybrid.cpp b/src/MOLECULE/angle_hybrid.cpp index 869b2bf41a..6cd641e73b 100644 --- a/src/MOLECULE/angle_hybrid.cpp +++ b/src/MOLECULE/angle_hybrid.cpp @@ -282,6 +282,16 @@ void AngleHybrid::coeff(int narg, char **arg) } } +/* ---------------------------------------------------------------------- + run angle style specific initialization +------------------------------------------------------------------------- */ + +void AngleHybrid::init_style() +{ + for (int m = 0; m < nstyles; m++) + if (styles[m]) styles[m]->init_style(); +} + /* ---------------------------------------------------------------------- return an equilbrium angle length ------------------------------------------------------------------------- */ diff --git a/src/MOLECULE/angle_hybrid.h b/src/MOLECULE/angle_hybrid.h index 9671e3a6cd..8f12bcb3db 100644 --- a/src/MOLECULE/angle_hybrid.h +++ b/src/MOLECULE/angle_hybrid.h @@ -27,11 +27,16 @@ namespace LAMMPS_NS { class AngleHybrid : public Angle { public: + int nstyles; // # of different angle styles + Angle **styles; // class list for each Angle style + char **keywords; // keyword for each Angle style + AngleHybrid(class LAMMPS *); ~AngleHybrid(); void compute(int, int); void settings(int, char **); void coeff(int, char **); + void init_style(); double equilibrium_angle(int); void write_restart(FILE *); void read_restart(FILE *); @@ -39,9 +44,6 @@ class AngleHybrid : public Angle { double memory_usage(); private: - int nstyles; // # of different angle styles - Angle **styles; // class list for each Angle style - char **keywords; // keyword for each Angle style int *map; // which style each angle type points to int *nanglelist; // # of angles in sub-style anglelists diff --git a/src/MOLECULE/dihedral_hybrid.h b/src/MOLECULE/dihedral_hybrid.h index 551a49d775..9c4b734996 100644 --- a/src/MOLECULE/dihedral_hybrid.h +++ b/src/MOLECULE/dihedral_hybrid.h @@ -27,6 +27,10 @@ namespace LAMMPS_NS { class DihedralHybrid : public Dihedral { public: + int nstyles; // # of different dihedral styles + Dihedral **styles; // class list for each Dihedral style + char **keywords; // keyword for each dihedral style + DihedralHybrid(class LAMMPS *); ~DihedralHybrid(); void compute(int, int); @@ -38,9 +42,6 @@ class DihedralHybrid : public Dihedral { double memory_usage(); private: - int nstyles; // # of different dihedral styles - Dihedral **styles; // class list for each Dihedral style - char **keywords; // keyword for each dihedral style int *map; // which style each dihedral type points to int *ndihedrallist; // # of dihedrals in sub-style dihedrallists diff --git a/src/MOLECULE/improper_hybrid.h b/src/MOLECULE/improper_hybrid.h index ceb75d0100..8e866fbf21 100644 --- a/src/MOLECULE/improper_hybrid.h +++ b/src/MOLECULE/improper_hybrid.h @@ -27,6 +27,10 @@ namespace LAMMPS_NS { class ImproperHybrid : public Improper { public: + int nstyles; // # of different improper styles + Improper **styles; // class list for each Improper style + char **keywords; // keyword for each improper style + ImproperHybrid(class LAMMPS *); ~ImproperHybrid(); void compute(int, int); @@ -37,9 +41,6 @@ class ImproperHybrid : public Improper { double memory_usage(); private: - int nstyles; // # of different improper styles - Improper **styles; // class list for each Improper style - char **keywords; // keyword for each improper style int *map; // which style each improper type points to int *nimproperlist; // # of impropers in sub-style improperlists diff --git a/src/USER-EWALDN/ewald_n.cpp b/src/USER-EWALDN/ewald_n.cpp index ed05ca5044..f3d7fe08d3 100644 --- a/src/USER-EWALDN/ewald_n.cpp +++ b/src/USER-EWALDN/ewald_n.cpp @@ -87,7 +87,6 @@ void EwaldN::init() error->all(FLERR,"Incorrect boundaries with slab EwaldN"); } - qqrd2e = force->qqrd2e; // check pair_style scale = 1.0; //mumurd2e = force->mumurd2e; //dielectric = force->dielectric; @@ -368,9 +367,11 @@ void EwaldN::init_self() memset(energy_self, 0, EWALD_NFUNCS*sizeof(double)); // self energy memset(virial_self, 0, EWALD_NFUNCS*sizeof(double)); + const double qscale = force->qqrd2e * scale; + if (function[0]) { // 1/r - virial_self[0] = -0.5*M_PI*qqrd2e*scale/(g2*volume)*sum[0].x*sum[0].x; - energy_self[0] = sum[0].x2*qqrd2e*scale*g1/sqrt(M_PI)-virial_self[0]; + virial_self[0] = -0.5*M_PI*qscale/(g2*volume)*sum[0].x*sum[0].x; + energy_self[0] = sum[0].x2*qscale*g1/sqrt(M_PI)-virial_self[0]; } if (function[1]) { // geometric 1/r^6 virial_self[1] = M_PI*sqrt(M_PI)*g3/(6.0*volume)*sum[1].x*sum[1].x; @@ -484,8 +485,9 @@ void EwaldN::compute_force() complex *cek, zc, zx, zxy; double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL; double *mu = atom->mu ? atom->mu[0] : NULL; + const double qscale = force->qqrd2e * scale; double *ke, c[EWALD_NFUNCS] = { - 8.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume), + 8.0*M_PI*qscale/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume), 2.0*M_PI*sqrt(M_PI)/(192.0*volume), 8.0*M_PI*mumurd2e/volume}; double kt = 4.0*pow(g_ewald, 3.0)/3.0/sqrt(M_PI)/c[3]; int i, kx, ky, lbytes = (2*nbox+1)*sizeof(cvector), *type = atom->type; @@ -587,8 +589,9 @@ void EwaldN::compute_energy(int eflag) complex *cek = cek_global; double *ke = kenergy; + const double qscale = force->qqrd2e * scale; double c[EWALD_NFUNCS] = { - 4.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume), + 4.0*M_PI*qscale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume), 2.0*M_PI*sqrt(M_PI)/(192.0*volume), 4.0*M_PI*mumurd2e/volume}; double sum[EWALD_NFUNCS]; int func[EWALD_NFUNCS]; @@ -625,8 +628,9 @@ void EwaldN::compute_virial(int vflag) complex *cek = cek_global; double *kv = kvirial; + const double qscale = force->qqrd2e * scale; double c[EWALD_NFUNCS] = { - 4.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume), + 4.0*M_PI*qscale/volume, 2.0*M_PI*sqrt(M_PI)/(24.0*volume), 2.0*M_PI*sqrt(M_PI)/(192.0*volume), 4.0*M_PI*mumurd2e/volume}; shape sum[EWALD_NFUNCS]; int func[EWALD_NFUNCS]; @@ -685,14 +689,16 @@ void EwaldN::compute_slabcorr(int eflag) while ((x+=3)qqrd2e * scale; - double ffact = -4.0*M_PI*qqrd2e*scale*dipole_all/volume; + double ffact = -4.0*M_PI*qscale*dipole_all/volume; double *f = atom->f[0]-1, *fn = f+3*atom->nlocal-3; q = atom->q; while ((f+=3)