int nbond_smoothcut( REAL_T *,  REAL_T *,  REAL_T *,  REAL_T *, REAL_T, REAL_T );

int nbond_smoothcut( REAL_T *x,  REAL_T *f,  REAL_T *enb,  REAL_T *eel, REAL_T enbfac, REAL_T eelfac )
/*
REAL_T *x, *f;
REAL_T *enb, *eel;
REAL_T enbfac, eelfac;
*/
{
	int iaci, iexcl, jexcl, jexcl_last, nexcl, skip;
	int ic, ll;
	REAL_T qi, qiqj;
	int i, j;
	REAL_T xi, yi, zi;
	REAL_T xij, yij, zij;
	REAL_T r2, rinv, r2inv, r6inv, f6, f12, dedx, dedy, dedz;

	REAL_T eelij;
	REAL_T de;
	REAL_T dij, erest;
	REAL_T cut2, icut2, rc, icut, icut6, icut12, rc_H;
	REAL_T A, B;

	*eel = *enb = 0.0;
	iexcl = 0;
	erest = 0.0;

	/* various powers of the cutoff calculated beforehand */
	cut2 = cut * cut;
	icut = 1. / cut;
	icut2 = 1. / cut2;
	icut6 = icut2 * icut2 * icut2;
	icut12 = icut6 * icut6;

#define CX 1.125

	/*what follows is mostly borrowed from egb */

	for (i = 0; i < prm->Natom; i++) {
		xi = x[3 * i];
		yi = x[3 * i + 1];
		zi = x[3 * i + 2];
		qi = prm->Charges[i];
		iaci = prm->Ntypes * (prm->Iac[i] - 1);
		jexcl = iexcl;
		jexcl_last = iexcl + prm->Iblo[i] - 1;
		nexcl = prm->ExclAt[jexcl] - 1;
		if (jexcl > jexcl_last)
			nexcl = -1;

		for (j = i + 1; j < prm->Natom; j++) {
			if (frozen[i] && frozen[j])
				continue;

			skip = 0;
			if (j == nexcl) {
				skip = 1;
				jexcl++;
				nexcl = prm->ExclAt[jexcl] - 1;
				if (jexcl > jexcl_last)
					nexcl = -1;
			}
			de = 0.0;
			if (!skip) {
				xij = xi - x[3 * j];
				yij = yi - x[3 * j + 1];
				zij = zi - x[3 * j + 2];
				qiqj = qi * prm->Charges[j];

				r2 = xij * xij + yij * yij + zij * zij;

				if (r2 < cut2) {

					dij = sqrt(r2);
					rinv = 1. / dij;
					r2inv = rinv * rinv;
					rc = (r2 * icut2) * (r2 * icut2) * (r2 * icut2);
					rc_H = rc * rc * (r2 * icut2) * (r2 * icut2);

					eelij =
						qiqj * (rinv + rc_H * CX * icut -
								rc_H * dij * icut2 - CX * icut);
					*eel += eelij;
					de -=
						qiqj * (rinv * r2inv - 18. * rc_H * r2inv * icut +
								17. * rc_H * rinv * icut2);
					/* smoothly get E_coul and derivative to 0 at cutoff */
					ic = prm->Cno[iaci + prm->Iac[j] - 1] - 1;

					if (ic > -1) {
						r6inv = r2inv * r2inv * r2inv;
						A = prm->Cn1[ic];
						B = prm->Cn2[ic];
						f6 = B * (r6inv - icut6);
						f12 = A * (r6inv * r6inv - icut12);
						*enb +=
							f12 - f6 + (2. * A * icut12 - B * icut6) * rc -
							3. * A * icut12 + 2. * B * icut6;
						de -=
							(12. * f12 - 6. * f6 - 12. * A * rc * icut12 +
							 6. * B * rc * icut6) * r2inv;
					}
					dedx = de * xij;
					dedy = de * yij;
					dedz = de * zij;
					f[3 * i] += dedx;
					f[3 * i + 1] += dedy;
					f[3 * i + 2] += dedz;
					f[3 * j] -= dedx;
					f[3 * j + 1] -= dedy;
					f[3 * j + 2] -= dedz;

				}				/* cutoff-if closed */
			}
		}

		iexcl += prm->Iblo[i];

	}
	return(0);
}