set up analytical long-range coulomb to have double precision accuracy (2826449b) · Commits · 郑智淋 / lammps

src/KSPACE/pair_lj_charmm_coul_long.cpp

+14 −21

Original line number	Diff line number	Diff line
		@@ -30,18 +30,14 @@
		#include "neighbor.h"
		#include "neigh_list.h"
		#include "neigh_request.h"
		#include "math_special.h"
		#include "math_const.h"
		#include "memory.h"
		#include "error.h"

		using namespace LAMMPS_NS;

		#define EWALD_F 1.12837917
		#define EWALD_P 0.3275911
		#define A1 0.254829592
		#define A2 -0.284496736
		#define A3 1.421413741
		#define A4 -1.453152027
		#define A5 1.061405429
		using namespace MathSpecial;
		using namespace MathConst;

		/* ---------------------------------------------------------------------- */

		@@ -144,11 +140,10 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
		if (!ncoultablebits \|\| rsq <= tabinnersq) {
		r = sqrt(rsq);
		grij = g_ewald * r;
		expm2 = exp(-grij*grij);
		t = 1.0 / (1.0 + EWALD_P*grij);
		erfc = t * (A1+t(A2+t(A3+t(A4+tA5)))) * expm2;
		expm2 = expmsq(grij);
		erfc = my_erfcx(grij) * expm2;
		prefactor = qqrd2e * qtmp*q[j]/r;
		forcecoul = prefactor * (erfc + EWALD_Fgrijexpm2);
		forcecoul = prefactor * (erfc + MY_ISPI4grijexpm2);
		if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
		} else {
		union_int_float_t rsq_lookup;
		@@ -453,11 +448,10 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
		if (!ncoultablebits \|\| rsq <= tabinnersq) {
		r = sqrt(rsq);
		grij = g_ewald * r;
		expm2 = exp(-grij*grij);
		t = 1.0 / (1.0 + EWALD_P*grij);
		erfc = t * (A1+t(A2+t(A3+t(A4+tA5)))) * expm2;
		expm2 = expmsq(grij);
		erfc = my_erfcx(grij) * expm2;
		prefactor = qqrd2e * qtmp*q[j]/r;
		forcecoul = prefactor * (erfc + EWALD_Fgrijexpm2 - 1.0);
		forcecoul = prefactor * (erfc + MY_ISPI4grijexpm2 - 1.0);
		if (rsq > cut_in_off_sq) {
		if (rsq < cut_in_on_sq) {
		rsw = (r - cut_in_off)/cut_in_diff;
		@@ -546,7 +540,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
		if (vflag) {
		if (rsq < cut_coulsq) {
		if (!ncoultablebits \|\| rsq <= tabinnersq) {
		forcecoul = prefactor * (erfc + EWALD_Fgrijexpm2);
		forcecoul = prefactor * (erfc + MY_ISPI4grijexpm2);
		if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
		} else {
		table = vtable[itable] + fraction*dvtable[itable];
		@@ -954,11 +948,10 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
		if (!ncoultablebits \|\| rsq <= tabinnersq) {
		r = sqrt(rsq);
		grij = g_ewald * r;
		expm2 = exp(-grij*grij);
		t = 1.0 / (1.0 + EWALD_P*grij);
		erfc = t * (A1+t(A2+t(A3+t(A4+tA5)))) * expm2;
		expm2 = expmsq(grij);
		erfc = my_erfcx(grij) * expm2;
		prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
		forcecoul = prefactor * (erfc + EWALD_Fgrijexpm2);
		forcecoul = prefactor * (erfc + MY_ISPI4grijexpm2);
		if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
		} else {
		union_int_float_t rsq_lookup;

src/KSPACE/pppm.cpp

+18 −14

Original line number	Diff line number	Diff line
		@@ -1209,23 +1209,24 @@ double PPPM::compute_qopt()
		sum2 = 0.0;
		sum3 = 0.0;
		sum4 = 0.0;
		const double inv_gew = 1.0/g_ewald;
		for (nx = -2; nx <= 2; nx++) {
		qx = unitkx(kper+nx_pppmnx);
		sx = exp(-0.25*square(qx/g_ewald));
		sx = expmsq(0.5qxinv_gew);
		argx = 0.5qxxprd/nx_pppm;
		wx = powsinxx(argx,twoorder);
		qx *= qx;

		for (ny = -2; ny <= 2; ny++) {
		qy = unitky(lper+ny_pppmny);
		sy = exp(-0.25*square(qy/g_ewald));
		sy = expmsq(0.5qyinv_gew);
		argy = 0.5qyyprd/ny_pppm;
		wy = powsinxx(argy,twoorder);
		qy *= qy;

		for (nz = -2; nz <= 2; nz++) {
		qz = unitkz(mper+nz_pppmnz);
		sz = exp(-0.25*square(qz/g_ewald));
		sz = expmsq(0.5qzinv_gew);
		argz = 0.5qzzprd_slab/nz_pppm;
		wz = powsinxx(argz,twoorder);
		qz *= qz;
		@@ -1297,7 +1298,7 @@ double PPPM::newton_raphson_f()
		double zprd = domain->zprd;
		bigint natoms = atom->natoms;

		double df_rspace = 2.0q2exp(-g_ewaldg_ewaldcutoff*cutoff) /
		double df_rspace = 2.0q2expmsq(g_ewald*cutoff) /
		sqrt(natomscutoffxprdyprdzprd);

		double df_kspace = compute_df_kspace();
		@@ -1338,7 +1339,7 @@ double PPPM::final_accuracy()

		double df_kspace = compute_df_kspace();
		double q2_over_sqrt = q2 / sqrt(natomscutoffxprdyprdzprd);
		double df_rspace = 2.0 * q2_over_sqrt * exp(-g_ewaldg_ewaldcutoff*cutoff);
		double df_rspace = 2.0 * q2_over_sqrt * expmsq(g_ewald*cutoff);
		double df_table = estimate_table_accuracy(q2_over_sqrt,df_rspace);
		double estimated_accuracy = sqrt(df_kspacedf_kspace + df_rspacedf_rspace +
		df_table*df_table);
		@@ -1579,22 +1580,23 @@ void PPPM::compute_gf_ik()
		numerator = 12.5663706/sqk;
		denominator = gf_denom(snx,sny,snz);
		sum1 = 0.0;
		const double inv_gew = 1.0/g_ewald;

		for (nx = -nbx; nx <= nbx; nx++) {
		qx = unitkx(kper+nx_pppmnx);
		sx = exp(-0.25*square(qx/g_ewald));
		sx = expmsq(0.5qxinv_gew);
		argx = 0.5qxxprd/nx_pppm;
		wx = powsinxx(argx,twoorder);

		for (ny = -nby; ny <= nby; ny++) {
		qy = unitky(lper+ny_pppmny);
		sy = exp(-0.25*square(qy/g_ewald));
		sy = expmsq(0.5qyinv_gew);
		argy = 0.5qyyprd/ny_pppm;
		wy = powsinxx(argy,twoorder);

		for (nz = -nbz; nz <= nbz; nz++) {
		qz = unitkz(mper+nz_pppmnz);
		sz = exp(-0.25*square(qz/g_ewald));
		sz = expmsq(0.5qzinv_gew);
		argz = 0.5qzzprd_slab/nz_pppm;
		wz = powsinxx(argz,twoorder);

		@@ -1662,6 +1664,7 @@ void PPPM::compute_gf_ik_triclinic()
		numerator = 12.5663706/sqk;
		denominator = gf_denom(snx,sny,snz);
		sum1 = 0.0;
		const double inv_gew = 1.0/g_ewald;

		for (nx = -nbx; nx <= nbx; nx++) {
		argx = MY_PIkper/nx_pppm + MY_PInx;
		@@ -1682,13 +1685,13 @@ void PPPM::compute_gf_ik_triclinic()
		x2lamdaT(&b[0],&b[0]);

		qx = unitk_lamda[0]+b[0];
		sx = exp(-0.25*square(qx/g_ewald));
		sx = expmsq(0.5qxinv_gew);

		qy = unitk_lamda[1]+b[1];
		sy = exp(-0.25*square(qy/g_ewald));
		sy = expmsq(0.5qyinv_gew);

		qz = unitk_lamda[2]+b[2];
		sz = exp(-0.25*square(qz/g_ewald));
		sz = expmsq(0.5qzinv_gew);

		dot1 = unitk_lamda[0]qx + unitk_lamda[1]qy + unitk_lamda[2]*qz;
		dot2 = qxqx+qyqy+qz*qz;
		@@ -1725,6 +1728,7 @@ void PPPM::compute_gf_ad()
		int k,l,m,n,kper,lper,mper;

		const int twoorder = 2*order;
		const double inv_gew = 1.0/g_ewald;

		for (int i = 0; i < 6; i++) sf_coeff[i] = 0.0;

		@@ -1733,7 +1737,7 @@ void PPPM::compute_gf_ad()
		mper = m - nz_pppm(2m/nz_pppm);
		qz = unitkz*mper;
		snz = square(sin(0.5qzzprd_slab/nz_pppm));
		sz = exp(-0.25*square(qz/g_ewald));
		sz = expmsq(0.5qzinv_gew);
		argz = 0.5qzzprd_slab/nz_pppm;
		wz = powsinxx(argz,twoorder);

		@@ -1741,7 +1745,7 @@ void PPPM::compute_gf_ad()
		lper = l - ny_pppm(2l/ny_pppm);
		qy = unitky*lper;
		sny = square(sin(0.5qyyprd/ny_pppm));
		sy = exp(-0.25*square(qy/g_ewald));
		sy = expmsq(0.5qyinv_gew);
		argy = 0.5qyyprd/ny_pppm;
		wy = powsinxx(argy,twoorder);

		@@ -1749,7 +1753,7 @@ void PPPM::compute_gf_ad()
		kper = k - nx_pppm(2k/nx_pppm);
		qx = unitkx*kper;
		snx = square(sin(0.5qxxprd/nx_pppm));
		sx = exp(-0.25*square(qx/g_ewald));
		sx = expmsq(0.5qxinv_gew);
		argx = 0.5qxxprd/nx_pppm;
		wx = powsinxx(argx,twoorder);

Admin message