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

Merge pull request #3482 from stanmoore1/slabcorr

Browse Source 
Fix issue with KSpace slab correction energy with non-neutral systems
This commit is contained in:
Axel Kohlmeyer
2022-10-11 21:12:38 -04:00
committed by GitHub
parent d499c2f61b 1e6321a1e9
commit 77740c4f07
20 changed files with 50 additions and 50 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/DIELECTRIC/pppm_dielectric.cpp
View File

@ -412,7 +412,7 @@ void PPPMDielectric::slabcorr()
double **x = atom->x;
double *eps = atom->epsilon;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
@ -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

6
src/DIELECTRIC/pppm_disp_dielectric.cpp
View File

@ -754,7 +754,7 @@ void PPPMDispDielectric::slabcorr(int /*eflag*/)
double *q = atom->q;
double **x = atom->x;
double *eps = atom->epsilon;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
@ -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

2
src/ELECTRODE/boundary_correction.h
View File

@ -27,7 +27,7 @@ class BoundaryCorrection : protected Pointers {
BoundaryCorrection(LAMMPS *);
virtual void vector_corr(double *, int, int, bool){};
virtual void matrix_corr(bigint *, double **){};
virtual void compute_corr(double, int, int, double &, double *){};
virtual void compute_corr(double, double, int, int, double &, double *){};
void setup(double, double, double);
void setup(double, double, double, double);

2
src/ELECTRODE/ewald_electrode.cpp
View File

@ -418,7 +418,7 @@ void EwaldElectrode::compute(int eflag, int vflag)
for (int j = 0; j < 6; j++) vatom[i][j] *= q[i] * qscale;
}
boundcorr->compute_corr(qsum, eflag_atom, eflag_global, energy, eatom);
boundcorr->compute_corr(qsum, slab_volfactor, eflag_atom, eflag_global, energy, eatom);
}
/* ---------------------------------------------------------------------- */

2
src/ELECTRODE/pppm_electrode.cpp
View File

@ -552,7 +552,7 @@ void PPPMElectrode::compute(int eflag, int vflag)
}
}
boundcorr->compute_corr(qsum, eflag_atom, eflag_global, energy, eatom);
boundcorr->compute_corr(qsum, slab_volfactor, eflag_atom, eflag_global, energy, eatom);
compute_vector_called = false;
}

4
src/ELECTRODE/slab_2d.cpp
View File

@ -32,8 +32,8 @@ using namespace MathConst;
------------------------------------------------------------------------- */
Slab2d::Slab2d(LAMMPS *lmp) : BoundaryCorrection(lmp){};
void Slab2d::compute_corr(double /*qsum*/, int eflag_atom, int eflag_global, double &energy,
double *eatom)
void Slab2d::compute_corr(double /*qsum*/, double /*slab_volfactor*/, int eflag_atom,
int eflag_global, double &energy, double *eatom)
{
double *q = atom->q;
double **x = atom->x;

2
src/ELECTRODE/slab_2d.h
View File

@ -27,7 +27,7 @@ class Slab2d : public BoundaryCorrection {
Slab2d(LAMMPS *);
void vector_corr(double *, int, int, bool) override;
void matrix_corr(bigint *, double **) override;
void compute_corr(double, int, int, double &, double *) override;
void compute_corr(double, double, int, int, double &, double *) override;
void setup(double);
};

10
src/ELECTRODE/slab_dipole.cpp
View File

@ -37,13 +37,13 @@ using namespace MathConst;
*/
SlabDipole::SlabDipole(LAMMPS *lmp) : BoundaryCorrection(lmp){};
void SlabDipole::compute_corr(double qsum, int eflag_atom, int eflag_global, double &energy,
double *eatom)
void SlabDipole::compute_corr(double qsum, double slab_volfactor, int eflag_atom,
int eflag_global, double &energy, double *eatom)
{
// compute local contribution to global dipole moment
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
for (int i = 0; i < nlocal; i++) dipole += q[i] * x[i][2];
@ -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

2
src/ELECTRODE/slab_dipole.h
View File

@ -27,7 +27,7 @@ class SlabDipole : public BoundaryCorrection {
SlabDipole(LAMMPS *);
void vector_corr(double *, int, int, bool);
void matrix_corr(bigint *, double **);
void compute_corr(double, int, int, double &, double *);
void compute_corr(double, double, int, int, double &, double *);
void setup(double);
};

4
src/ELECTRODE/wire_dipole.cpp
View File

@ -34,8 +34,8 @@ using namespace MathConst;
*/
WireDipole::WireDipole(LAMMPS *lmp) : BoundaryCorrection(lmp){};
void WireDipole::compute_corr(double /*qsum*/, int eflag_atom, int eflag_global, double &energy,
double *eatom)
void WireDipole::compute_corr(double /*qsum*/, double /*slab_volfactor*/, int eflag_atom,
int eflag_global, double &energy, double *eatom)
{
double *q = atom->q;
double **x = atom->x;

2
src/ELECTRODE/wire_dipole.h
View File

@ -27,7 +27,7 @@ class WireDipole : public BoundaryCorrection {
WireDipole(LAMMPS *);
void vector_corr(double *, int, int, bool);
void matrix_corr(bigint *, double **);
void compute_corr(double, int, int, double &, double *);
void compute_corr(double, double, int, int, double &, double *);
void setup(double);
};

2
src/INTEL/pppm_electrode_intel.cpp
View File

@ -272,7 +272,7 @@ void PPPMElectrodeIntel::compute(int eflag, int vflag)
slabflag = 0; // bypass compute_second's slabcorr()
PPPMIntel::compute_second(eflag, vflag);
slabflag = tempslabflag;
boundcorr->compute_corr(qsum, eflag_atom, eflag_global, energy, eatom);
boundcorr->compute_corr(qsum, slab_volfactor, eflag_atom, eflag_global, energy, eatom);
compute_vector_called = false;
}

6
src/KOKKOS/pppm_kokkos.cpp
View File

@ -2760,7 +2760,7 @@ void PPPMKokkos<DeviceType>::slabcorr()
{
// compute local contribution to global dipole moment
zprd = domain->zprd;
zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atomKK->nlocal;
double dipole = 0.0;
@ -2791,7 +2791,7 @@ void PPPMKokkos<DeviceType>::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<DeviceType>::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<class DeviceType>

10
src/KSPACE/ewald.cpp
View File

@ -1171,7 +1171,7 @@ void Ewald::slabcorr()
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
@ -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
@ -1380,7 +1380,7 @@ void Ewald::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag)
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -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

6
src/KSPACE/ewald_disp.cpp
View File

@ -1426,7 +1426,7 @@ void EwaldDisp::compute_slabcorr()
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
@ -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

10
src/KSPACE/pppm.cpp
View File

@ -2915,7 +2915,7 @@ void PPPM::slabcorr()
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
@ -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
@ -3431,7 +3431,7 @@ void PPPM::slabcorr_groups(int groupbit_A, int groupbit_B, int AA_flag)
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int *mask = atom->mask;
int nlocal = atom->nlocal;
@ -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

6
src/KSPACE/pppm_cg.cpp
View File

@ -595,7 +595,7 @@ void PPPMCG::slabcorr()
const double * const q = atom->q;
const double * const * const x = atom->x;
const double zprd = domain->zprd;
const double zprd_slab = domain->zprd*slab_volfactor;
double dipole = 0.0;
@ -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);
}
}

6
src/KSPACE/pppm_disp.cpp
View File

@ -8110,7 +8110,7 @@ void PPPMDisp::slabcorr(int /*eflag*/)
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
double dipole = 0.0;
@ -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

6
src/KSPACE/pppm_disp_tip4p.cpp
View File

@ -498,7 +498,7 @@ void PPPMDispTIP4P::slabcorr(int /*eflag*/)
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
int *type = atom->type;
double *xi, xM[3]; int iH1, iH2; //for TIP4P virtual site
@ -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

6
src/KSPACE/pppm_tip4p.cpp
View File

@ -490,7 +490,7 @@ void PPPMTIP4P::slabcorr()
double *q = atom->q;
double **x = atom->x;
double zprd = domain->zprd;
double zprd_slab = domain->zprd*slab_volfactor;
int nlocal = atom->nlocal;
int *type = atom->type;
double *xi, xM[3]; int iH1, iH2; //for TIP4P virtual site
@ -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