Fix bug in ewald_dipole forces (b9e33e63) · Commits · 郑智淋 / lammps

src/KSPACE/ewald_dipole.cpp

+41 −43

Original line number	Diff line number	Diff line
		@@ -443,7 +443,7 @@ void EwaldDipole::compute(int eflag, int vflag)
		sypz = sn[ky][1][i]cs[kz][2][i] + cs[ky][1][i]sn[kz][2][i];
		exprl = cs[kx][0][i]cypz - sn[kx][0][i]sypz;
		expim = sn[kx][0][i]cypz + cs[kx][0][i]sypz;
		partial = expimsfacrl_all[k] - exprlsfacim_all[k];
		partial = (muk[k][i])(expimsfacrl_all[k] - exprl*sfacim_all[k]);
		ek[i][0] += partial*eg[k][0];
		ek[i][1] += partial*eg[k][1];
		ek[i][2] += partial*eg[k][2];
		@@ -465,14 +465,12 @@ void EwaldDipole::compute(int eflag, int vflag)
		const double muscale = qqrd2e * scale;

		for (i = 0; i < nlocal; i++) {
		for (k = 0; k < kcount; k++) {
		//f[i][0] += qscale * q[i]*ek[i][0];
		//f[i][1] += qscale * q[i]*ek[i][1];
		//if (slabflag != 2) f[i][2] += qscale * q[i]*ek[i][2];
		f[i][0] += muscale * muk[k][i] * ek[i][0];
		f[i][1] += muscale * muk[k][i] * ek[i][1];
		if (slabflag != 2) f[i][2] += muscale * muk[k][i] * ek[i][2];
		}
		f[i][0] += muscale * ek[i][0];
		f[i][1] += muscale * ek[i][1];
		if (slabflag != 2) f[i][2] += muscale * ek[i][2];
		}

		// sum global energy across Kspace vevs and add in volume-dependent term
		@@ -548,7 +546,7 @@ void EwaldDipole::eik_dot_r()
		// define first val. of cos and sin

		for (ic = 0; ic < 3; ic++) {
		sqk = (unitk[ic] * unitk[ic]);
		sqk = unitk[ic]*unitk[ic];
		if (sqk <= gsqmx) {
		cstr1 = 0.0;
		sstr1 = 0.0;
		@@ -557,8 +555,8 @@ void EwaldDipole::eik_dot_r()
		sn[0][ic][i] = 0.0;
		cs[1][ic][i] = cos(unitk[ic]*x[i][ic]);
		sn[1][ic][i] = sin(unitk[ic]*x[i][ic]);
		cs[-1][ic][i] = cs[1][0][i];
		sn[-1][ic][i] = -sn[1][0][i];
		cs[-1][ic][i] = cs[1][ic][i];
		sn[-1][ic][i] = -sn[1][ic][i];
		muk[n][i] = (mu[i][ic]*unitk[ic]);
		cstr1 += muk[n][i]*cs[1][ic][i];
		sstr1 += muk[n][i]*sn[1][ic][i];

Admin message