ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@11756 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2014-04-09 16:50:23 +00:00
parent 5e3ff6cffd
commit c0fbdbefaa
14 changed files with 107 additions and 170 deletions
Download Patch File Download Diff File Expand all files Collapse all files

46
src/KSPACE/ewald.cpp
View File

@ -106,8 +106,7 @@ void Ewald::init()
// extract short-range Coulombic cutoff from pair style
scale = 1.0;
triclinic = domain->triclinic;
pair_check();
int itmp;
@ -116,25 +115,12 @@ void Ewald::init()
error->all(FLERR,"KSpace style is incompatible with Pair style");
double cutoff = *p_cutoff;
qsum = qsqsum = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
qsum += atom->q[i];
qsqsum += atom->q[i]*atom->q[i];
}
// compute qsum & qsqsum and warn if not charge-neutral
double tmp;
MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum = tmp;
MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsqsum = tmp;
if (qsqsum == 0.0)
error->all(FLERR,"Cannot use kspace solver on system with no charge");
if (fabs(qsum) > SMALL && comm->me == 0) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->warning(FLERR,str);
}
scale = 1.0;
qqrd2e = force->qqrd2e;
qsum_qsq(0);
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -143,11 +129,8 @@ void Ewald::init()
// setup K-space resolution
q2 = qsqsum * force->qqrd2e;
bigint natoms = atom->natoms;
triclinic = domain->triclinic;
// use xprd,yprd,zprd even if triclinic so grid size is the same
// adjust z dimension for 2d slab Ewald
// 3d Ewald just uses zprd since slab_volfactor = 1.0
@ -440,7 +423,7 @@ void Ewald::compute(int eflag, int vflag)
// convert E-field to force
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
for (i = 0; i < nlocal; i++) {
f[i][0] += qscale * q[i]*ek[i][0];
@ -448,12 +431,19 @@ void Ewald::compute(int eflag, int vflag)
if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2];
}
// global energy
// sum global energy across Kspace vevs and add in volume-dependent term
// reset qsum and qsqsum if atom count has changed
if (eflag_global) {
for (k = 0; k < kcount; k++)
energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
sfacim_all[k]*sfacim_all[k]);
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qscale;
@ -1206,7 +1196,7 @@ void Ewald::slabcorr()
const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) energy += qscale * e_slabcorr;
@ -1337,7 +1327,7 @@ void Ewald::compute_group_group(int groupbit_A, int groupbit_B, int AA_flag)
MPI_Allreduce(sfacrl_B,sfacrl_B_all,kcount,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(sfacim_B,sfacim_B_all,kcount,MPI_DOUBLE,MPI_SUM,world);
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
double partial_group;
// total group A <--> group B energy
@ -1437,7 +1427,7 @@ void Ewald::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag)
// compute corrections
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
const double efact = qscale * MY_2PI/volume;
e2group += efact * (dipole_A*dipole_B - 0.5*(qsum_A*dipole_r2_B +

51
src/KSPACE/msm.cpp
View File

@ -166,46 +166,23 @@ void MSM::init()
error->all(FLERR,"Cannot (yet) use single precision with MSM "
"(remove -DFFT_SINGLE from Makefile and recompile)");
pair_check();
// extract short-range Coulombic cutoff from pair style
triclinic = domain->triclinic;
pair_check();
int itmp;
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
if (p_cutoff == NULL)
error->all(FLERR,"KSpace style is incompatible with Pair style");
cutoff = *p_cutoff;
// compute qsum & qsqsum and give error if not charge-neutral
// compute qsum & qsqsum and error if not charge-neutral
qsum = qsqsum = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
qsum += atom->q[i];
qsqsum += atom->q[i]*atom->q[i];
}
qqrd2e = force->qqrd2e;
scale = 1.0;
double tmp;
MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum = tmp;
MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsqsum = tmp;
q2 = qsqsum * force->qqrd2e;
if (qsqsum == 0.0)
error->all(FLERR,"Cannot use kspace solver on system with no charge");
// not yet sure of the correction needed for non-neutral systems
if (fabs(qsum) > SMALL) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->all(FLERR,str);
}
triclinic = domain->triclinic;
qqrd2e = force->qqrd2e;
qsum_qsq(1);
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -586,16 +563,22 @@ void MSM::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// total long-range energy
// sum global energy across procs and add in self-energy term
// reset qsum and qsqsum if atom count has changed
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
double e_self = qsqsum*gamma(0.0)/cutoff; // Self-energy term
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
double e_self = qsqsum*gamma(0.0)/cutoff;
energy -= e_self;
energy *= 0.5*qscale;
}
@ -2796,7 +2779,7 @@ void MSM::fieldforce()
// convert E-field to force
const double qfactor = force->qqrd2e*scale*q[i];
const double qfactor = qqrd2e*scale*q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
f[i][2] += qfactor*ekz;

1
src/KSPACE/msm.h
View File

@ -40,7 +40,6 @@ class MSM : public KSpace {
double precision;
int nfactors;
int *factors;
double qsum,qsqsum,q2;
double qqrd2e;
double cutoff;
double volume;

46
src/KSPACE/pppm.cpp
View File

@ -48,8 +48,8 @@ using namespace MathSpecial;
#define MAXORDER 7
#define OFFSET 16384
#define SMALL 0.00001
#define LARGE 10000.0
#define SMALL 0.00001
#define EPS_HOC 1.0e-7
enum{REVERSE_RHO};
@ -203,8 +203,6 @@ void PPPM::init()
// extract short-range Coulombic cutoff from pair style
triclinic = domain->triclinic;
scale = 1.0;
pair_check();
int itmp = 0;
@ -250,26 +248,10 @@ void PPPM::init()
// compute qsum & qsqsum and warn if not charge-neutral
qsum = qsqsum = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
qsum += atom->q[i];
qsqsum += atom->q[i]*atom->q[i];
}
double tmp;
MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum = tmp;
MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsqsum = tmp;
q2 = qsqsum * force->qqrd2e;
if (qsqsum == 0.0)
error->all(FLERR,"Cannot use kspace solver on system with no charge");
if (fabs(qsum) > SMALL && me == 0) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->warning(FLERR,str);
}
scale = 1.0;
qqrd2e = force->qqrd2e;
qsum_qsq(0);
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -678,14 +660,20 @@ void PPPM::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// sum global energy across procs and add in volume-dependent term
// reset qsum and qsqsum if atom count has changed
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@ -2460,7 +2448,7 @@ void PPPM::fieldforce_ik()
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
const double qfactor = qqrd2e * scale * q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
if (slabflag != 2) f[i][2] += qfactor*ekz;
@ -2531,7 +2519,7 @@ void PPPM::fieldforce_ad()
// convert E-field to force and substract self forces
const double qfactor = force->qqrd2e * scale;
const double qfactor = qqrd2e * scale;
s1 = x[i][0]*hx_inv;
s2 = x[i][1]*hy_inv;
@ -2954,7 +2942,7 @@ void PPPM::slabcorr()
const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) energy += qscale * e_slabcorr;
@ -3141,7 +3129,7 @@ void PPPM::compute_group_group(int groupbit_A, int groupbit_B, int AA_flag)
poisson_groups(AA_flag);
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
// total group A <--> group B energy
// self and boundary correction terms are in compute_group_group.cpp
@ -3487,7 +3475,7 @@ void PPPM::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag)
// compute corrections
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
const double efact = qscale * MY_2PI/volume;
e2group += efact * (dipole_A*dipole_B - 0.5*(qsum_A*dipole_r2_B +

1
src/KSPACE/pppm.h
View File

@ -53,7 +53,6 @@ class PPPM : public KSpace {
int me,nprocs;
int nfactors;
int *factors;
double qsum,qsqsum,q2;
double cutoff;
double volume;
double delxinv,delyinv,delzinv,delvolinv;

14
src/KSPACE/pppm_cg.cpp
View File

@ -209,14 +209,20 @@ void PPPMCG::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// sum global energy across procs and add in volume-dependent term
// reset qsum and qsqsum if atom count has changed
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@ -406,7 +412,7 @@ void PPPMCG::fieldforce_ik()
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
const double qfactor = qqrd2e * scale * q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
if (slabflag != 2) f[i][2] += qfactor*ekz;
@ -479,7 +485,7 @@ void PPPMCG::fieldforce_ad()
// convert E-field to force and substract self forces
const double qfactor = force->qqrd2e * scale;
const double qfactor = qqrd2e * scale;
s1 = x[i][0]*hx_inv;
s2 = x[i][1]*hy_inv;
@ -618,7 +624,7 @@ void PPPMCG::slabcorr()
const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) energy += qscale * e_slabcorr;

53
src/KSPACE/pppm_disp.cpp
View File

@ -233,14 +233,9 @@ void PPPMDisp::init()
deallocate();
deallocate_peratom();
// set scale
scale = 1.0;
triclinic = domain->triclinic;
// check whether cutoff and pair style are set
triclinic = domain->triclinic;
pair_check();
int tmp;
@ -269,7 +264,8 @@ void PPPMDisp::init()
case 1:
k = 0; break;
case 6:
if ((ewald_mix==GEOMETRIC || ewald_mix==SIXTHPOWER|| mixflag == 1) && mixflag!= 2) { k = 1; break; }
if ((ewald_mix==GEOMETRIC || ewald_mix==SIXTHPOWER ||
mixflag == 1) && mixflag!= 2) { k = 1; break; }
else if (ewald_mix==ARITHMETIC && mixflag!=2) { k = 2; break; }
else if (mixflag == 2) { k = 3; break; }
default:
@ -282,41 +278,20 @@ void PPPMDisp::init()
// warn, if function[0] is not set but charge attribute is set!
if (!function[0] && atom->q_flag && me == 0) {
char str[128];
sprintf(str, "Charges are set, but coulombic solver is not used");
error->warning(FLERR, str);
}
// compute qsum & qsqsum, if function[0] is set, print error if no charges are set or warn if not charge-neutral
// compute qsum & qsqsum, if function[0] is set, warn if not charge-neutral
scale = 1.0;
qqrd2e = force->qqrd2e;
natoms_original = atom->natoms;
if (function[0]) {
if (!atom->q_flag)
error->all(FLERR,"Kspace style with selected options "
"requires atom attribute q");
qsum = qsqsum = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
qsum += atom->q[i];
qsqsum += atom->q[i]*atom->q[i];
}
double tmp;
MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum = tmp;
MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsqsum = tmp;
if (qsqsum == 0.0)
error->all(FLERR,"Cannot use kspace solver with selected options "
"on system with no charge");
if (fabs(qsum) > SMALL && me == 0) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->warning(FLERR,str);
}
}
if (function[0]) qsum_qsq(0);
// if kspace is TIP4P, extract TIP4P params from pair style
// bond/angle are not yet init(), so insure equilibrium request is valid
@ -351,8 +326,8 @@ void PPPMDisp::init()
alpha = qdist / (cos(0.5*theta) * blen);
}
// initialize the pair style to get the coefficients
neighrequest_flag = 0;
pair->init();
neighrequest_flag = 1;
@ -1144,6 +1119,7 @@ void PPPMDisp::compute(int eflag, int vflag)
}
// sum energy across procs and add in volume-dependent term
// reset qsum and qsqsum if atom count has changed
const double qscale = force->qqrd2e * scale;
if (eflag_global) {
@ -1156,6 +1132,11 @@ void PPPMDisp::compute(int eflag, int vflag)
energy_1 *= 0.5*volume;
energy_6 *= 0.5*volume;
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
energy_1 -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy_6 += - MY_PI*MY_PIS/(6*volume)*pow(g_ewald_6,3)*csumij +

1
src/KSPACE/pppm_disp.h
View File

@ -62,7 +62,6 @@ Variables needed for calculating the 1/r and 1/r^6 potential
int me,nprocs;
int nfactors;
int *factors;
double qsum,qsqsum;
double csumij;
double csum;
double *csumi; //needed as correction term for per atom calculations!

1
src/KSPACE/pppm_disp_tip4p.cpp
View File

@ -97,7 +97,6 @@ void PPPMDispTIP4P::particle_map_c(double delx, double dely, double delz,
p2g[i][1] = ny;
p2g[i][2] = nz;
// check that entire stencil around nx,ny,nz will fit in my 3d brick
if (nx+nlow < nxlo || nx+nup > nxhi ||

41
src/KSPACE/pppm_old.cpp
View File

@ -154,8 +154,6 @@ void PPPMOld::init()
// extract short-range Coulombic cutoff from pair style
scale = 1.0;
pair_check();
int itmp=0;
@ -198,25 +196,10 @@ void PPPMOld::init()
// compute qsum & qsqsum and warn if not charge-neutral
qsum = qsqsum = 0.0;
for (int i = 0; i < atom->nlocal; i++) {
qsum += atom->q[i];
qsqsum += atom->q[i]*atom->q[i];
}
double tmp;
MPI_Allreduce(&qsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsum = tmp;
MPI_Allreduce(&qsqsum,&tmp,1,MPI_DOUBLE,MPI_SUM,world);
qsqsum = tmp;
if (qsqsum == 0.0)
error->all(FLERR,"Cannot use kspace solver on system with no charge");
if (fabs(qsum) > SMALL && me == 0) {
char str[128];
sprintf(str,"System is not charge neutral, net charge = %g",qsum);
error->warning(FLERR,str);
}
scale = 1.0;
qqrd2e = force->qqrd2e;
qsum_qsq(0);
natoms_original = atom->natoms;
// set accuracy (force units) from accuracy_relative or accuracy_absolute
@ -742,14 +725,20 @@ void PPPMOld::compute(int eflag, int vflag)
if (evflag_atom) fieldforce_peratom();
// sum global energy across procs and add in volume-dependent term
// reset qsum and qsqsum if atom count has changed
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@ -971,8 +960,6 @@ void PPPMOld::set_grid()
acons[7][5] = 1755948832039.0 / 36229939200000.0;
acons[7][6] = 4887769399.0 / 37838389248.0;
double q2 = qsqsum * force->qqrd2e;
// use xprd,yprd,zprd even if triclinic so grid size is the same
// adjust z dimension for 2d slab PPPM
// 3d PPPM just uses zprd since slab_volfactor = 1.0
@ -2232,7 +2219,7 @@ void PPPMOld::fieldforce()
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
const double qfactor = qqrd2e * scale * q[i];
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
if (slabflag != 2) f[i][2] += qfactor*ekz;
@ -2471,7 +2458,7 @@ void PPPMOld::slabcorr()
const double e_slabcorr = MY_2PI*(dipole_all*dipole_all -
qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume;
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) energy += qscale * e_slabcorr;
@ -2645,7 +2632,7 @@ void PPPMOld::compute_group_group(int groupbit_A, int groupbit_B, int BA_flag)
poisson_groups(BA_flag);
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
// total group A <--> group B energy
// self and boundary correction terms are in compute_group_group.cpp

1
src/KSPACE/pppm_old.h
View File

@ -52,7 +52,6 @@ class PPPMOld : public KSpace {
int me,nprocs;
int nfactors;
int *factors;
double qsum,qsqsum;
double cutoff;
double volume;
double delxinv,delyinv,delzinv,delvolinv;

13
src/KSPACE/pppm_stagger.cpp
View File

@ -203,14 +203,20 @@ void PPPMStagger::compute(int eflag, int vflag)
}
// sum global energy across procs and add in volume-dependent term
// reset qsum and qsqsum if atom count has changed
const double qscale = force->qqrd2e * scale;
const double qscale = qqrd2e * scale;
if (eflag_global) {
double energy_all;
MPI_Allreduce(&energy,&energy_all,1,MPI_DOUBLE,MPI_SUM,world);
energy = energy_all;
if (atom->natoms != natoms_original) {
qsum_qsq(0);
natoms_original = atom->natoms;
}
energy *= 0.5*volume/float(nstagger);
energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
@ -813,7 +819,7 @@ void PPPMStagger::fieldforce_ik()
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i] / float(nstagger);
const double qfactor = qqrd2e * scale * q[i] / float(nstagger);
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
if (slabflag != 2) f[i][2] += qfactor*ekz;
@ -884,7 +890,7 @@ void PPPMStagger::fieldforce_ad()
// convert E-field to force and substract self forces
const double qfactor = force->qqrd2e * scale / float(nstagger);
const double qfactor = qqrd2e * scale / float(nstagger);
s1 = x[i][0]*hx_inv + stagger;
s2 = x[i][1]*hy_inv + stagger;
@ -899,7 +905,6 @@ void PPPMStagger::fieldforce_ad()
sf *= 2*q[i]*q[i];
f[i][1] += qfactor*(eky*q[i] - sf);
sf = sf_coeff[4]*sin(2*MY_PI*s3);
sf += sf_coeff[5]*sin(4*MY_PI*s3);
sf *= 2*q[i]*q[i];

3
src/KSPACE/pppm_stagger.h
View File

@ -52,7 +52,8 @@ class PPPMStagger : public PPPM {
inline double gf_denom2(const double &x, const double &y,
const double &z) const {
const double &z) const
{
double sx,sy,sz;
double x2 = x*x;
double y2 = y*y;

5
src/KSPACE/pppm_tip4p.cpp
View File

@ -223,7 +223,7 @@ void PPPMTIP4P::fieldforce_ik()
// convert E-field to force
const double qfactor = force->qqrd2e * scale * q[i];
const double qfactor = qqrd2e * scale * q[i];
if (type[i] != typeO) {
f[i][0] += qfactor*ekx;
f[i][1] += qfactor*eky;
@ -326,7 +326,8 @@ void PPPMTIP4P::fieldforce_ad()
ekz *= hz_inv;
// convert E-field to force and substract self forces
const double qfactor = force->qqrd2e * scale;
const double qfactor = qqrd2e * scale;
s1 = x[i][0]*hx_inv;
s2 = x[i][1]*hy_inv;