From 877f4e7e29cc7624fcc8118ff95e8a09e4e6cbae Mon Sep 17 00:00:00 2001 From: Stan Gerald Moore Date: Tue, 11 Oct 2022 09:43:18 -0600 Subject: [PATCH] Fix issue with KSpace slab correction energy with non-neutral systems --- src/DIELECTRIC/pppm_dielectric.cpp | 4 ++-- src/DIELECTRIC/pppm_disp_dielectric.cpp | 4 ++-- src/ELECTRODE/slab_dipole.cpp | 4 ++-- src/KOKKOS/pppm_kokkos.cpp | 4 ++-- src/KSPACE/ewald.cpp | 6 +++--- src/KSPACE/ewald_disp.cpp | 4 ++-- src/KSPACE/pppm.cpp | 6 +++--- src/KSPACE/pppm_cg.cpp | 4 ++-- src/KSPACE/pppm_disp.cpp | 4 ++-- src/KSPACE/pppm_disp_tip4p.cpp | 4 ++-- src/KSPACE/pppm_tip4p.cpp | 4 ++-- 11 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/DIELECTRIC/pppm_dielectric.cpp b/src/DIELECTRIC/pppm_dielectric.cpp index c063750a75..cf2b19008c 100644 --- a/src/DIELECTRIC/pppm_dielectric.cpp +++ b/src/DIELECTRIC/pppm_dielectric.cpp @@ -442,7 +442,7 @@ void PPPMDielectric::slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -453,7 +453,7 @@ void PPPMDielectric::slabcorr() double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * eps[i]*q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp index beb392a3b3..e07e8dcb01 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.cpp +++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp @@ -793,7 +793,7 @@ void PPPMDispDielectric::slabcorr(int /*eflag*/) // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -804,7 +804,7 @@ void PPPMDispDielectric::slabcorr(int /*eflag*/) double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * eps[i]*q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections diff --git a/src/ELECTRODE/slab_dipole.cpp b/src/ELECTRODE/slab_dipole.cpp index 2e04a012fc..0e49f5c134 100644 --- a/src/ELECTRODE/slab_dipole.cpp +++ b/src/ELECTRODE/slab_dipole.cpp @@ -65,7 +65,7 @@ void SlabDipole::compute_corr(double qsum, int eflag_atom, int eflag_global, dou // compute corrections double const e_slabcorr = MY_2PI * - (dipole_all * dipole_all - qsum * dipole_r2 - qsum * qsum * zprd * zprd / 12.0) / volume; + (dipole_all * dipole_all - qsum * dipole_r2 - qsum * qsum * zprd_slab * zprd_slab / 12.0) / volume; double const qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -75,7 +75,7 @@ void SlabDipole::compute_corr(double qsum, int eflag_atom, int eflag_global, dou for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i] * (x[i][2] * dipole_all - 0.5 * (dipole_r2 + qsum * x[i][2] * x[i][2]) - - qsum * zprd * zprd / 12.0); + qsum * zprd_slab * zprd_slab / 12.0); } // add on force corrections diff --git a/src/KOKKOS/pppm_kokkos.cpp b/src/KOKKOS/pppm_kokkos.cpp index 432e866b2c..5171c398ad 100644 --- a/src/KOKKOS/pppm_kokkos.cpp +++ b/src/KOKKOS/pppm_kokkos.cpp @@ -2791,7 +2791,7 @@ void PPPMKokkos::slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -2833,7 +2833,7 @@ KOKKOS_INLINE_FUNCTION void PPPMKokkos::operator()(TagPPPM_slabcorr3, const int &i) const { d_eatom[i] += efact * q[i]*(x(i,2)*dipole_all - 0.5*(dipole_r2 + - qsum*x(i,2)*x(i,2)) - qsum*zprd*zprd/12.0); + qsum*x(i,2)*x(i,2)) - qsum*zprd_slab*zprd_slab/12.0); } template diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp index 2f787c54fe..f0638a16e0 100644 --- a/src/KSPACE/ewald.cpp +++ b/src/KSPACE/ewald.cpp @@ -1200,7 +1200,7 @@ void Ewald::slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -1211,7 +1211,7 @@ void Ewald::slabcorr() double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections @@ -1436,7 +1436,7 @@ void Ewald::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag) const double efact = qscale * MY_2PI/volume; e2group += efact * (dipole_A*dipole_B - 0.5*(qsum_A*dipole_r2_B + - qsum_B*dipole_r2_A) - qsum_A*qsum_B*zprd*zprd/12.0); + qsum_B*dipole_r2_A) - qsum_A*qsum_B*zprd_slab*zprd_slab/12.0); // add on force corrections diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp index d5ebc4dd36..4d1f27fb65 100644 --- a/src/KSPACE/ewald_disp.cpp +++ b/src/KSPACE/ewald_disp.cpp @@ -1465,7 +1465,7 @@ void EwaldDisp::compute_slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = force->qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -1476,7 +1476,7 @@ void EwaldDisp::compute_slabcorr() double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index e94193759f..7dcecd8238 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -2944,7 +2944,7 @@ void PPPM::slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -2955,7 +2955,7 @@ void PPPM::slabcorr() double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections @@ -3487,7 +3487,7 @@ void PPPM::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag) const double efact = qscale * MY_2PI/volume; e2group += efact * (dipole_A*dipole_B - 0.5*(qsum_A*dipole_r2_B + - qsum_B*dipole_r2_A) - qsum_A*qsum_B*zprd*zprd/12.0); + qsum_B*dipole_r2_A) - qsum_A*qsum_B*zprd_slab*zprd_slab/12.0); // add on force corrections diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp index 077eb9f3f4..44518ed056 100644 --- a/src/KSPACE/pppm_cg.cpp +++ b/src/KSPACE/pppm_cg.cpp @@ -629,7 +629,7 @@ void PPPMCG::slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = qqrd2e * scale; if (eflag_global) energy += qscale * e_slabcorr; @@ -641,7 +641,7 @@ void PPPMCG::slabcorr() for (j = 0; j < num_charged; j++) { i = is_charged[j]; eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } } diff --git a/src/KSPACE/pppm_disp.cpp b/src/KSPACE/pppm_disp.cpp index 7377ec6f39..cf62ef4520 100644 --- a/src/KSPACE/pppm_disp.cpp +++ b/src/KSPACE/pppm_disp.cpp @@ -8139,7 +8139,7 @@ void PPPMDisp::slabcorr(int /*eflag*/) // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = force->qqrd2e * scale; if (eflag_global) energy_1 += qscale * e_slabcorr; @@ -8150,7 +8150,7 @@ void PPPMDisp::slabcorr(int /*eflag*/) double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections diff --git a/src/KSPACE/pppm_disp_tip4p.cpp b/src/KSPACE/pppm_disp_tip4p.cpp index 9ebc32ad05..76f3fabfc3 100644 --- a/src/KSPACE/pppm_disp_tip4p.cpp +++ b/src/KSPACE/pppm_disp_tip4p.cpp @@ -534,7 +534,7 @@ void PPPMDispTIP4P::slabcorr(int /*eflag*/) // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = force->qqrd2e * scale; if (eflag_global) energy_1 += qscale * e_slabcorr; @@ -545,7 +545,7 @@ void PPPMDispTIP4P::slabcorr(int /*eflag*/) double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp index e9d72bcefe..76ca14a30e 100644 --- a/src/KSPACE/pppm_tip4p.cpp +++ b/src/KSPACE/pppm_tip4p.cpp @@ -526,7 +526,7 @@ void PPPMTIP4P::slabcorr() // compute corrections const double e_slabcorr = MY_2PI*(dipole_all*dipole_all - - qsum*dipole_r2 - qsum*qsum*zprd*zprd/12.0)/volume; + qsum*dipole_r2 - qsum*qsum*zprd_slab*zprd_slab/12.0)/volume; const double qscale = force->qqrd2e * scale; if (eflag_global) energy_1 += qscale * e_slabcorr; @@ -537,7 +537,7 @@ void PPPMTIP4P::slabcorr() double efact = qscale * MY_2PI/volume; for (int i = 0; i < nlocal; i++) eatom[i] += efact * q[i]*(x[i][2]*dipole_all - 0.5*(dipole_r2 + - qsum*x[i][2]*x[i][2]) - qsum*zprd*zprd/12.0); + qsum*x[i][2]*x[i][2]) - qsum*zprd_slab*zprd_slab/12.0); } // add on force corrections