[Code scan] Clear pppm/dplr electric forces when charges vanish

[njzjz]

Found during a Codex global scan of deepmodeling/deepmd-kit at commit 73de44b1f94471b2e3bdb6b11f57b34d7bc791bb.

Problem

pppm/dplr can leave stale electric-field forces in fele when all charges become zero.

Evidence:

• PPPMDPLR::init() allocates and clears fele once:

  deepmd-kit/source/lmp/pppm_dplr.cpp

  Lines 51 to 64 in 73de44b

  	void PPPMDPLR::init() {
  	// DPLR PPPM requires newton on, b/c it computes forces on ghost atoms
  	if (force->newton == 0) {
  	error->all(FLERR, "Kspace style pppm/dplr requires newton on");
  	}
  	PPPM::init();
  	int nlocal = atom->nlocal;
  	// cout << " ninit pppm/dplr ---------------------- " << nlocal << endl;
  	fele.resize(static_cast<size_t>(nlocal) * 3);
  	fill(fele.begin(), fele.end(), 0.0);
  	}

• compute() returns early when qsqsum == 0.0:

  deepmd-kit/source/lmp/pppm_dplr.cpp

  Lines 82 to 93 in 73de44b

  	// if atom count has changed, update qsum and qsqsum
  	if (atom->natoms != natoms_original) {
  	qsum_qsq();
  	natoms_original = atom->natoms;
  	}
  	// return if there are no charges
  	if (qsqsum == 0.0) {
  	return;
  	}

• fele is cleared only inside fieldforce_ik() and fieldforce_ad(), which are not reached on the early return:

  deepmd-kit/source/lmp/pppm_dplr.cpp

  Lines 280 to 302 in 73de44b

  	void PPPMDPLR::fieldforce_ik() {
  	int i, l, m, n, nx, ny, nz, mx, my, mz;
  	FFT_SCALAR dx, dy, dz, x0, y0, z0;
  	FFT_SCALAR ekx, eky, ekz;
  	// loop over my charges, interpolate electric field from nearby grid points
  	// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
  	// (dx,dy,dz) = distance to "lower left" grid pt
  	// (mx,my,mz) = global coords of moving stencil pt
  	// ek = 3 components of E-field on particle
  	double* q = atom->q;
  	double** x = atom->x;
  	// double **f = atom->f;
  	int nlocal = atom->nlocal;
  	int nghost = atom->nghost;
  	int nall = nlocal + nghost;
  	fele.resize(static_cast<size_t>(nlocal) * 3);
  	fill(fele.begin(), fele.end(), 0.0);
  	for (i = 0; i < nlocal; i++) {

  and

  deepmd-kit/source/lmp/pppm_dplr.cpp

  Lines 344 to 377 in 73de44b

  	void PPPMDPLR::fieldforce_ad() {
  	int i, l, m, n, nx, ny, nz, mx, my, mz;
  	FFT_SCALAR dx, dy, dz;
  	FFT_SCALAR ekx, eky, ekz;
  	double s1, s2, s3;
  	double sf = 0.0;
  	double* prd;
  	prd = domain->prd;
  	double xprd = prd[0];
  	double yprd = prd[1];
  	double zprd = prd[2];
  	double hx_inv = nx_pppm / xprd;
  	double hy_inv = ny_pppm / yprd;
  	double hz_inv = nz_pppm / zprd;
  	// loop over my charges, interpolate electric field from nearby grid points
  	// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
  	// (dx,dy,dz) = distance to "lower left" grid pt
  	// (mx,my,mz) = global coords of moving stencil pt
  	// ek = 3 components of E-field on particle
  	double* q = atom->q;
  	double** x = atom->x;
  	// double **f = atom->f;
  	int nlocal = atom->nlocal;
  	int nghost = atom->nghost;
  	int nall = nlocal + nghost;
  	fele.resize(static_cast<size_t>(nlocal) * 3);
  	fill(fele.begin(), fele.end(), 0.0);

• FixDPLR::post_force() reads pppm_dplr->get_fele() and adds it to the DPLR force input:

  deepmd-kit/source/lmp/fix_dplr.cpp

  Lines 637 to 644 in 73de44b

  	// revise force according to efield
  	if (pppm_dplr) {
  	const vector<double>& dfele_(pppm_dplr->get_fele());
  	assert(dfele_.size() == nlocal * 3);
  	for (int ii = 0; ii < nlocal * 3; ++ii) {
  	dfele[ii] += dfele_[ii];
  	}
  	}

Impact

If a simulation changes from nonzero charges to all-zero charges, pppm/dplr can expose the previous step's electric-field forces to fix dplr, producing a spurious DPLR correction after the Coulomb field should be zero.

Suggested Fix

Resize and clear fele before the qsqsum == 0.0 early return, or clear it in that branch. Add a charge-changing regression where the first step has nonzero charges and a later step has all charges zero, asserting fix dplr sees zero PPPM electric force.