/***=======================================================================***/
/*** Like the nonbonded loop, there are different versions of the bond,    ***/
/*** angle, and dihedral routines.  The different variants can calculate   ***/
/*** forces alone, energies alone, forces and energies, or forces,         ***/
/*** energies, and virial contributions.                                   ***/
/***=======================================================================***/
#if NEEDFORCE == 1
  #if NEEDENERGY == 1
    #if NEEDVIRIAL == 1 
      #define ATTNPFRC AttenuatePairVir
      #define BONDCALC BondVir
      #define ANGLCALC AnglVir
      #define DIHECALC DiheVir
    #else
      #define ATTNPFRC AttenuatePairFrcNrg
      #define BONDCALC BondFrcNrg
      #define ANGLCALC AnglFrcNrg
      #define DIHECALC DiheFrcNrg
    #endif
  #else
    #define ATTNPFRC AttenuatePairFrc
    #define BONDCALC BondFrc
    #define ANGLCALC AnglFrc
    #define DIHECALC DiheFrc
  #endif
#else
  #define ATTNPFRC AttenuatePairNrg
  #define BONDCALC BondNrg
  #define ANGLCALC AnglNrg
  #define DIHECALC DiheNrg
#endif

/***=======================================================================***/
/*** AttenuatePairForce: routine to modify the pair force between atoms in ***/
/***                     order to account for exclusion of basic van-der   ***/
/***                     Waals or electrostatic interactions.  After any   ***/
/***                     necessary corrections are added, the forces,      ***/
/***                     energies, and virials are accumulated.            ***/
/***                                                                       ***/
/*** Arguments:                                                            ***/
/***=======================================================================***/
#if NEEDFORCE == 1
#if NEEDENERGY == 1
void ATTNPFRC(prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc, int atmA, int atmB,
	      double dx, double dy, double dz, double fmag, double* afrc,
	      double* bfrc, Energy *sysUV, double elec14fac, double lj14fac)
#else
void ATTNPFRC(prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc, int atmA, int atmB,
	      double dx, double dy, double dz, double fmag, double* afrc,
	      double* bfrc, double elec14fac, double lj14fac)
#endif
#else
void ATTNPFRC(prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc, int atmA, int atmB,
	      double dx, double dy, double dz, Energy *sysUV, double elec14fac,
	      double lj14fac)
#endif
{
  int irc, irh;
  double r, r2, invr, invr2, qq;
#if NEEDFORCE == 1
  CSpln *CdSD, *HdSD;
#endif
#if NEEDENERGY == 1
  CSpln *CSD, *HSD;
#endif

  /*** (Re)compute r and related quantities ***/
  r2 = dx*dx + dy*dy + dz*dz;
  r = sqrt(r2);
  invr = 1.0/r;
  invr2 = invr*invr;

  /*** Correct 1:2 electrostatics.  All electrostatic interactions are ***/
  /*** counted in the nonbonded routine (direct and reciprocal space), ***/
  /*** so exclusions or 1:4 screened interactions must be addressed.   ***/
  /*** However, electrostatic direct-space interactions are computed   ***/
  /*** with a coarse lookup table, which is inaccurate at short range. ***/
  /*** If an electrostatic interaction was calculated at short range,  ***/
  /*** we subtract off the contributions from the coarse lookup table  ***/
  /*** and replace them with calculations from a more accurate fine    ***/
  /*** lookup table before subtracting the analytic qq/r interaction.  ***/
  qq = tp->Charges[atmA] * tp->Charges[atmB];
  if (r2 < MINNB2) {
    irc = r2*Cfrc->ivdr;
#if NEEDFORCE == 1
    CdSD = Cfrc->dSD;
    fmag -= qq*(((CdSD[irc].A*r2 + CdSD[irc].B)*r2 + CdSD[irc].C)*r2 +
		CdSD[irc].D);
#endif
#if NEEDENERGY == 1
    CSD = Cfrc->SD;
    sysUV->delec -= qq*(((CSD[irc].A*r2 + CSD[irc].B)*r2 + CSD[irc].C)*r2 +
			CSD[irc].D);
#endif
    irh = r2*Hfrc->ivdr;
#if NEEDFORCE == 1
    HdSD = Hfrc->dSD;
    fmag += qq*(((HdSD[irh].A*r2 + HdSD[irh].B)*r2 + HdSD[irh].C)*r2 +
		HdSD[irh].D);
#endif
#if NEEDENERGY == 1
    HSD = Hfrc->SD;
    sysUV->delec += qq*(((HSD[irh].A*r2 + HSD[irh].B)*r2 + HSD[irh].C)*r2 +
			HSD[irh].D);
#endif
  }
#if NEEDFORCE == 1
  fmag += elec14fac*BIOQ*qq*invr*invr2;
#endif
#if NEEDENERGY == 1
  sysUV->delec -= elec14fac*BIOQ*qq*invr;
#endif

  /*** Correct 1:2 van-der Waals interactions ***/
  if (r2 > MINNB2) {

    int ljA, ljB;
    double invr4, invr6;

    ljA = tp->LJIdx[atmA];
    ljB = tp->LJIdx[atmB];
    if (ljA >= 0 && ljB >= 0) {
      invr4 = invr2*invr2;
      invr6 = invr4*invr2;
#if NEEDFORCE == 1
      fmag -= lj14fac*invr4*invr4*(tp->LJftab.map[ljA][2*ljB]*invr6 +
				   tp->LJftab.map[ljA][2*ljB+1]);
#endif
#if NEEDENERGY == 1
      sysUV->vdw12 -= lj14fac*invr6*invr6*tp->LJutab.map[ljA][2*ljB];
      sysUV->vdw6 -= lj14fac*invr6*tp->LJutab.map[ljA][2*ljB+1];
#endif
    }
  }

#if NEEDFORCE == 1
  /*** Accumulate the force ***/
  afrc[0] += fmag*dx;
  afrc[1] += fmag*dy;
  afrc[2] += fmag*dz;
  bfrc[0] -= fmag*dx;
  bfrc[1] -= fmag*dy;
  bfrc[2] -= fmag*dz;

#if NEEDVIRIAL == 1
  /*** Accumulate the stress tensor ***/
  sysUV->Vir[0] += dx*fmag*dx;
  sysUV->Vir[1] += dx*fmag*dy;
  sysUV->Vir[2] += dx*fmag*dz;
  sysUV->Vir[3] += dy*fmag*dx;
  sysUV->Vir[4] += dy*fmag*dy;
  sysUV->Vir[5] += dy*fmag*dz;
  sysUV->Vir[6] += dz*fmag*dx;
  sysUV->Vir[7] += dz*fmag*dy;
  sysUV->Vir[8] += dz*fmag*dz;
#endif
#endif
}

/***=======================================================================***/
/*** BondFrc: computes the force on two atoms due to a bonded interaction. ***/
/***                                                                       ***/
/*** Arguments:                                                            ***/
/***   [ab]ptr:    the coordinates for atoms A or B                        ***/
/***   [ab]frc:    the forces for atoms A or B                             ***/
/***   bcom:       the bond command                                        ***/
/***   tp:         the topology                                            ***/
/***   Cfrc:       the "coarse" electrostatic force/energy spline table    ***/
/***   Hfrc:       the "fine" electrostatic force/energy spline table      ***/
/***   sysUV:      the system energy and virial (state information)        ***/
/***=======================================================================***/
#if NEEDFORCE == 1
#if NEEDENERGY == 1
void BONDCALC(double *aptr, double *bptr, double *afrc, double *bfrc,
	      bondcomm *bcom, prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc,
	      Energy *sysUV)
#else
void BONDCALC(double *aptr, double *bptr, double *afrc, double *bfrc,
              bondcomm *bcom, prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc)
#endif
#else
void BONDCALC(double *aptr, double *bptr, bondcomm *bcom, prmtop *tp,
              FrcTab *Cfrc, FrcTab *Hfrc, Energy *sysUV)
#endif
{
  double r, r2, dx, dy, dz, dl;
#if NEEDFORCE == 1
  double fmag;
#endif

  /*** Compute displacement ***/
  dx = bptr[0] - aptr[0];
  dy = bptr[1] - aptr[1];
  dz = bptr[2] - aptr[2];

  /*** Accumulate the bond force and energy ***/
  r2 = dx*dx + dy*dy + dz*dz;
  r = sqrt(r2);
  dl = tp->BParam[bcom->t].l0 - r;
#if NEEDFORCE == 1
  fmag = -2.0*tp->BParam[bcom->t].K*dl/r;
#endif
#if NEEDENERGY == 1
  sysUV->bond += tp->BParam[bcom->t].K*dl*dl;
  sysUV->BondUdc[3*bcom->t] += r*r;
  sysUV->BondUdc[3*bcom->t+1] += r*tp->BParam[bcom->t].l0;
  sysUV->BondUdc[3*bcom->t+2] += pow(tp->BParam[bcom->t].l0, 2.0);
#endif

#if NEEDFORCE == 1
#if NEEDVIRIAL == 0
#if NEEDENERGY == 1
  AttenuatePairFrcNrg(tp, Cfrc, Hfrc, bcom->a, bcom->b, dx, dy, dz, fmag,
		      afrc, bfrc, sysUV, 1.0, 1.0);
#else
  AttenuatePairFrc(tp, Cfrc, Hfrc, bcom->a, bcom->b, dx, dy, dz, fmag,
		   afrc, bfrc, 1.0, 1.0);
#endif
#else
  AttenuatePairVir(tp, Cfrc, Hfrc, bcom->a, bcom->b, dx, dy, dz, fmag,
                   afrc, bfrc, sysUV, 1.0, 1.0);
#endif
#else
  AttenuatePairNrg(tp, Cfrc, Hfrc, bcom->a, bcom->b, dx, dy, dz, sysUV, 1.0,
		   1.0);
#endif
}

/***=======================================================================***/
/*** AnglFrc: computes the force on two atoms due to an angle interaction. ***/
/***                                                                       ***/
/*** Arguments:                                                            ***/
/***   [abc]ptr:   the coordinates for atoms A, B, or C                    ***/
/***   [abc]frc:   the forces for atoms A, B, or C                         ***/
/***   acom:       the angle command                                       ***/
/***   tp:         the topology                                            ***/
/***   Cfrc:       the "coarse" electrostatic force/energy spline table    ***/
/***   Hfrc:       the "fine" electrostatic force/energy spline table      ***/
/***   sysUV:      the system energy and virial (state information)        ***/
/***=======================================================================***/
#if NEEDFORCE == 1
#if NEEDENERGY == 1
void ANGLCALC(double *aptr, double *bptr, double *cptr, double *afrc,
	      double *bfrc, double *cfrc, anglcomm *acom, prmtop *tp,
	      FrcTab *Cfrc, FrcTab *Hfrc, Energy *sysUV)
#else
void ANGLCALC(double *aptr, double *bptr, double *cptr, double *afrc,
	      double *bfrc, double *cfrc, anglcomm *acom, prmtop *tp,
	      FrcTab *Cfrc, FrcTab *Hfrc)
#endif
#else
void ANGLCALC(double *aptr, double *bptr, double *cptr, anglcomm *acom,
              prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc, Energy *sysUV)
#endif
{
  int i;
  double costheta, theta, mgba, mgbc;
#if NEEDFORCE == 1
  double dA, mbabc, sqba, sqbc, adf, cdf;
#endif
  double invbabc, dtheta;
  double ac[3], ba[3], bc[3];

  /*** Compute displacements ***/
  for (i = 0; i < 3; i++) {
    ba[i] = aptr[i] - bptr[i];
    bc[i] = cptr[i] - bptr[i];
    ac[i] = cptr[i] - aptr[i];
  }

  /*** On to the angle force computation ***/
  mgba = ba[0]*ba[0] + ba[1]*ba[1] + ba[2]*ba[2];
  mgbc = bc[0]*bc[0] + bc[1]*bc[1] + bc[2]*bc[2];
  invbabc = 1.0/sqrt(mgba*mgbc);
  costheta = (ba[0]*bc[0] + ba[1]*bc[1] + ba[2]*bc[2]) * invbabc;
  costheta = (costheta < -1.0) ? -1.0 : (costheta > 1.0) ? 1.0 : costheta;
  theta = acos(costheta);
  dtheta = theta - tp->AParam[acom->t].th0;
#if NEEDFORCE == 1
  dA = -2.0*tp->AParam[acom->t].K*dtheta /
    sqrt(1.0 - costheta*costheta);
  sqba = dA/mgba;
  sqbc = dA/mgbc;
  mbabc = dA * invbabc;

  /*** Accumulate the angle forces and stress tensor ***/
  for (i = 0; i < 3; i++) {
    adf = bc[i]*mbabc - costheta*ba[i]*sqba;
    cdf = ba[i]*mbabc - costheta*bc[i]*sqbc;
    afrc[i] -= adf;
    bfrc[i] += adf + cdf;
    cfrc[i] -= cdf;
#if NEEDVIRIAL == 1
    sysUV->Vir[0+i] += ba[0]*adf + bc[0]*cdf;
    sysUV->Vir[3+i] += ba[1]*adf + bc[1]*cdf;
    sysUV->Vir[6+i] += ba[2]*adf + bc[2]*cdf;
#endif
  }
#endif

#if NEEDENERGY == 1
  /*** Accumulate the angle energy ***/
  sysUV->angl += tp->AParam[acom->t].K*dtheta*dtheta;
#endif

  /*** Bail out if this angle's A:C interactions do not need to be ***/
  /*** excluded (i.e. they were already excluded by some bond)     ***/
  if (acom->excl == 0) {
    return; 
  }

#if NEEDFORCE == 1
#if NEEDVIRIAL == 0
#if NEEDENERGY == 1
  AttenuatePairFrcNrg(tp, Cfrc, Hfrc, acom->a, acom->c, ac[0], ac[1], ac[2],
		      0.0, afrc, cfrc, sysUV, 1.0, 1.0);
#else
  AttenuatePairFrc(tp, Cfrc, Hfrc, acom->a, acom->c, ac[0], ac[1], ac[2], 0.0,
                   afrc, cfrc, 1.0, 1.0);
#endif
#else
  AttenuatePairVir(tp, Cfrc, Hfrc, acom->a, acom->c, ac[0], ac[1], ac[2], 0.0,
                   afrc, cfrc, sysUV, 1.0, 1.0);
#endif
#else
  AttenuatePairNrg(tp, Cfrc, Hfrc, acom->a, acom->c, ac[0], ac[1], ac[2],
		   sysUV, 1.0, 1.0);
#endif

}

/***=======================================================================***/
/*** DiheFrc: computes the force on four atoms due to a dihedral motion.   ***/
/***          This also applies to improper dihedrals.  Standard and       ***/
/***          improper dihedrals are treated using the following notation: ***/
/***                                                                       ***/
/***             Standard Dihedral          Improper Dihedral              ***/
/***                                                                       ***/
/***                       D                           D                   ***/
/***                      /                           /                    ***/
/***                  B--C                        B--C                     ***/
/***                 /                                \                    ***/
/***                A                                  A                   ***/
/***                                                                       ***/
/*** Arguments:                                                            ***/
/***   [abcd]ptr:   the coordinates for atoms A, B, C, or D                ***/
/***   [abcd]frc:   the forces for atoms A, B, C, or D                     ***/
/***   hcom:        the dihedral command, giving details on all Fourier    ***/
/***                series terms and a flag to tell whether this is an     ***/
/***                improper or a standard dihedral                        ***/
/***   hdef:        the array of dihedral definitions                      ***/
/***=======================================================================***/
#if NEEDFORCE == 1
#if NEEDENERGY == 1
void DIHECALC(double *aptr, double *bptr, double *cptr, double *dptr,
	      double *afrc, double *bfrc, double *cfrc, double *dfrc,
	      dihecomm *hcom, prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc,
	      Energy *sysUV)
#else
void DIHECALC(double *aptr, double *bptr, double *cptr, double *dptr,
	      double *afrc, double *bfrc, double *cfrc, double *dfrc,
	      dihecomm *hcom, prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc)
#endif
#else
void DIHECALC(double *aptr, double *bptr, double *cptr, double *dptr,
	      dihecomm *hcom, prmtop *tp, FrcTab *Cfrc, FrcTab *Hfrc,
	      Energy *sysUV)
#endif
{
  int i, idx;
  double theta, costheta;
#if NEEDFORCE == 1
  double fr, fa, fb1, fc1, fb2, fc2, fd, isinb2, isinc2;
  double mgab, mgbc, mgcd, invab, invbc, invcd, invabc, invbcd, cosb, cosc;
#endif
  double sangle;
  double ab[3], bc[3], cd[3], crabbc[3], crbccd[3], scr[3];
  dihedef *hdef;

  /*** Compute displacements ***/
  for (i = 0; i < 3; i++) {
    ab[i] = bptr[i] - aptr[i];
    bc[i] = cptr[i] - bptr[i];
    cd[i] = dptr[i] - cptr[i];
  }
  CrossP(ab, bc, crabbc);
  CrossP(bc, cd, crbccd);
  costheta = crabbc[0]*crbccd[0] + crabbc[1]*crbccd[1] + crabbc[2]*crbccd[2];
  costheta /=
    sqrt((crabbc[0]*crabbc[0] + crabbc[1]*crabbc[1] + crabbc[2]*crabbc[2]) *
	 (crbccd[0]*crbccd[0] + crbccd[1]*crbccd[1] + crbccd[2]*crbccd[2]));
  CrossP(crabbc, crbccd, scr);
  costheta = (costheta < -1.0) ? -1.0 : (costheta > 1.0) ? 1.0 : costheta;
  if (scr[0]*bc[0] + scr[1]*bc[1] + scr[2]*bc[2] > 0.0) {
    theta = acos(costheta);
  }
  else {
    theta = -acos(costheta);
  }

  /*** Compute the magnitude of the force and accumulate the energy ***/
#if NEEDFORCE == 1
  fr = 0.0;
#endif
  hdef = tp->HParam;
  for (i = 0; i < hcom->nt; i++) {
    idx = hcom->t[i];
    sangle = hdef[idx].N*theta - hdef[idx].Phi;
#if NEEDFORCE == 1
    fr += hdef[idx].K * hdef[idx].N * sin(sangle);
#endif
#if NEEDENERGY == 1
    sysUV->dihe += hdef[idx].K * (1.0 + cos(sangle));
#endif
  }

#if NEEDFORCE == 1
  /*** Other pre-computations ***/
  mgab = sqrt(ab[0]*ab[0] + ab[1]*ab[1] + ab[2]*ab[2]);
  invab = 1.0/mgab;
  mgbc = sqrt(bc[0]*bc[0] + bc[1]*bc[1] + bc[2]*bc[2]);
  invbc = 1.0/mgbc;
  mgcd = sqrt(cd[0]*cd[0] + cd[1]*cd[1] + cd[2]*cd[2]);
  invcd = 1.0/mgcd;
  cosb = -(ab[0]*bc[0] + ab[1]*bc[1] + ab[2]*bc[2])*invab*invbc;
  isinb2 = (cosb*cosb < 0.9999) ? 1.0/(1.0 - cosb*cosb) : 0.0;
  cosc = -(bc[0]*cd[0] + bc[1]*cd[1] + bc[2]*cd[2])*invbc*invcd;
  isinc2 = (cosc*cosc < 0.9999) ? 1.0/(1.0 - cosc*cosc) : 0.0;
  isinb2 *= fr;
  isinc2 *= fr;
  invabc = invab*invbc;
  invbcd = invbc*invcd;
  for (i = 0; i < 3; i++) {
    crabbc[i] *= invabc;
    crbccd[i] *= invbcd;
  }

  /*** Transform the dihedral forces to cartesian coordinates ***/
  fa = -invab * isinb2;
  fb1 = (mgbc - mgab*cosb) * invabc * isinb2;
  fb2 = cosc * invbc * isinc2;
  fc1 = (mgbc - mgcd*cosc) * invbcd * isinc2;
  fc2 = cosb * invbc * isinb2;
  fd = -invcd * isinc2;
  for (i = 0; i < 3; i++) {
    afrc[i] += crabbc[i] * fa;
    bfrc[i] += fb1 * crabbc[i] - fb2 * crbccd[i];
    cfrc[i] += -fc1 * crbccd[i] + fc2 * crabbc[i];
    dfrc[i] += -fd * crbccd[i];
  }
#endif

  /*** Evaluate 1-4 interactions ***/
  if (hcom->eval14 == 1) {
#if NEEDFORCE == 1
#if NEEDVIRIAL == 0
#if NEEDENERGY == 1
  AttenuatePairFrcNrg(tp, Cfrc, Hfrc, hcom->a, hcom->d, dptr[0] - aptr[0],
		      dptr[1] - aptr[1], dptr[2] - aptr[2], 0.0, afrc, dfrc,
		      sysUV, hcom->scee, hcom->scnb);
#else
  AttenuatePairFrc(tp, Cfrc, Hfrc, hcom->a, hcom->d, dptr[0] - aptr[0],
		   dptr[1] - aptr[1], dptr[2] - aptr[2], 0.0, afrc, dfrc,
		   hcom->scee, hcom->scnb);
#endif
#else
  AttenuatePairVir(tp, Cfrc, Hfrc, hcom->a, hcom->d, dptr[0] - aptr[0],
		   dptr[1] - aptr[1], dptr[2] - aptr[2], 0.0, afrc, dfrc,
		   sysUV, hcom->scee, hcom->scnb);
#endif
#else
  AttenuatePairNrg(tp, Cfrc, Hfrc, hcom->a, hcom->d, dptr[0] - aptr[0],
		   dptr[1] - aptr[1], dptr[2] - aptr[2], sysUV, hcom->scee,
		   hcom->scnb);
#endif
  }
}

#undef BONDCALC
#undef ANGLCALC
#undef DIHECALC
#undef ATTNPFRC