/*
 * Calculate the generalized Born energy and first derivatives.
 *
 * Calling parameters are as follows:
 *
 * natom - input: the number of atoms
 * x - input: the atomic (x,y,z) coordinates
 * f - updated: the gradient vector
 * fs - input: overlap parameters
 * rborn - input: atomic radii
 * q - input: atomic charges
 * kappa - input: inverse of the Debye-Huckel length
 * diel_ext - input: solvent dielectric constant
 * enb - updated: Lennard-Jones energy
 * eelt - updated: gas-phase electrostatic energy
 */


REAL_T egb( int *natom, REAL_T *x, REAL_T *f, REAL_T *fs, REAL_T *rborn, REAL_T *q,
	    REAL_T *kappa, REAL_T *diel_ext, REAL_T *enb, REAL_T *eelt )

#define BOFFSET 0.09
#define KSCALE 0.73

#define PIx4 12.5663706143591724639918
#define PIx2  6.2831853071795862319959
#define PIx1  3.1415926535897931159979

{



  int i, i3, j, j3, k;
  int iexcl, jexcl, jexcl_last, nexcl, npairs, ic, iaci;
  int *iexw;
  REAL_T epol, dielfac, qi, qj, qiqj, fgbi, fgbk, rb2, expmkf;
  REAL_T elec, evdw, sumda, daix, daiy, daiz;
  REAL_T xi, yi, zi, xij, yij, zij;
  REAL_T dedx, dedy, dedz, de;
  REAL_T dij1i, dij3i, temp1;
  REAL_T qi2h, qid2h, datmp;
  REAL_T theta, ri1i, dij2i;

  REAL_T dij, sumi;
  REAL_T eel, f6, f12, rinv, r2inv, r6inv;
  REAL_T r2, ri, rj, sj, sj2, thi;
  REAL_T uij, efac, temp4, temp5, temp6;
  REAL_T *reff, *sumdeijda, *sumepol, *sumelec, *sumevdw, *psi;

  /* LCPO stuff follows */
  int count, count2, icount;
  REAL_T si, sumAij, sumAjk, sumAijAjk, sumdAijddijdxi;
  REAL_T sumdAijddijdyi, sumdAijddijdzi, sumdAijddijdxiAjk;
  REAL_T sumdAijddijdyiAjk, sumdAijddijdziAjk, rij, tmpaij, Aij,
    dAijddij;
  REAL_T dAijddijdxj, dAijddijdyj, dAijddijdzj;
  REAL_T sumdAjkddjkdxj, sumdAjkddjkdyj, sumdAjkddjkdzj, p3p4Aij;
  REAL_T xk, yk, zk, rjk2, djk1i, rjk, vdw2dif, tmpajk, Ajk, sumAjk2,
    dAjkddjk;
  REAL_T dAjkddjkdxj, dAjkddjkdyj, dAjkddjkdzj, lastxj, lastyj, lastzj;
  REAL_T dAidxj, dAidyj, dAidzj, Ai, dAidxi, dAidyi,
    dAidzi;
  REAL_T totsasa;
  REAL_T xj, yj, zj;
  REAL_T dumbo,tmpsd;
  REAL_T rgbmax1i, rgbmax2i, rgbmaxpsmax2;

  /* arrays and vars for use in vectorizable loops */
  int icount2;
  REAL_T *vfgbi, *vefac,*r2x,*vrinv,*vexpmkf,*vectmp1,*vectmp2,*vectmp3,*vectmp4;
  REAL_T * vectmp5, *vectmp6, *vectmp7;
  REAL_T kkscale,invdiel_ext;
  int *nb1,*nbl1,*nbl2,*nbl3,*nbl4,*nbl5;
  int k1,k2,k3,k4,k5;
  REAL_T *ivd1,*ivd2,*ivd3,*ivd4,*ivd5;
  REAL_T *vsj,*vsj1,*vsj2,*vsj3,*vsj4,*vsj5;
  REAL_T *d1,*d2,*vuij,*vlog1;
  REAL_T logrgbmax1i,sumitmp;
  

  /*FGB taylor coefficients follow */
  /* from A to H :                 */
  /* 1/3 , 2/5 , 3/7 , 4/9 , 5/11  */
  /* 4/3 , 12/5 , 24/7 , 40/9 , 60/11 */

#define TA 0.33333333333333333333
#define TB 0.4
#define TC 0.42857142857142857143
#define TD 0.44444444444444444444
#define TDD 0.45454545454545454545

#define TE 1.33333333333333333333
#define TF 2.4
#define TG 3.42857142857142857143
#define TH 4.44444444444444444444
#define THH 5.45454545454545454545

  /*
   * Smooth "cut-off" in calculating GB effective radii.
   * Implementd by Andreas Svrcek-Seiler and Alexey Onufriev.
   * The integration over solute is performed up to rgbmax and includes
   * parts of spheres; that is an atom is not just "in" or "out", as
   * with standard non-bonded cut.  As a result, calclated effective
   * radii are less than rgbmax. This saves time, and there is no
   * discontinuity in dReff/drij.
   *
   * Only the case rgbmax > 5*max(sij) = 5*fsmax ~ 9A is handled; this is
   * enforced in mdread().  Smaller values would not make much physical
   * sense anyway.
   *
   * Note: rgbmax must be less than or equal to cut so that the pairlist
   * generated from cut may be applied to calculation of the effective
   * radius and its derivatives.
   */

  if (rgbmax > cut) {
    printf("Error in egb: rgbmax = %f is greater than cutoff = %f\n",
	   rgbmax, cut);
    exit(1);
  }
  rgbmax1i = 1.0/rgbmax;
  rgbmax2i = rgbmax1i*rgbmax1i;
  rgbmaxpsmax2 = (rgbmax+prm->Fsmax)*(rgbmax+prm->Fsmax);
  
  logrgbmax1i = log(rgbmax1i);
  kkscale = -(*kappa)*KSCALE;
  invdiel_ext = 1./(*diel_ext);
 
  /* allocate arrays for use in vectorizable loops */

  vfgbi = vector(0, *natom);
  vefac = vector(0, *natom);
  r2x = vector(0, *natom);
  vrinv = vector(0, *natom);
  vexpmkf = vector(0, *natom);
  vectmp1 = vector(0, *natom);
  vectmp2 = vector(0, *natom);
  vectmp3 = vector(0, *natom);
  vectmp4 = vector(0, *natom);
  vectmp5 = vector(0, *natom);
  vectmp6 = vector(0, *natom);
  vectmp7 = vector(0, *natom);
   
  nb1 = ivector(0,*natom);
  nbl1 = ivector(0,*natom);
  nbl2 = ivector(0,*natom);
  nbl3 = ivector(0,150);
  nbl4 = ivector(0,50);
  nbl5 = ivector(0,5);
  
  
  ivd1 = vector(0, *natom);
  ivd2 = vector(0, *natom);
  ivd3 =  vector(0, 150);
  ivd4 =  vector(0, 50);
  ivd5 = vector(0, 5);

  vsj = vector(0, *natom);
  vsj1 = vector(0, *natom);
  vsj2 = vector(0, *natom);
  vsj3 =  vector(0, 150);
  vsj4 =  vector(0, 50);
  vsj5 = vector(0, 5);

  d1 = vector(0, *natom);
  d2 = vector(0, *natom);
  vuij = vector(0, *natom);
  vlog1 = vector(0, *natom);
  
  /*  
  REAL_T *ivd1,*ivd2,*ivd3,*ivd4,*ivd5;
  REAL_T *vsj1,*vsj2,*vsj3,*vsj4,*vsj5;
  REAL_T *d1,*d2,*vuij,*vlog1;
  */
  sumepol = vector(0, *natom);
  sumelec = vector(0, *natom);
  sumevdw = vector(0, *natom);

  /* Allocate and initialize the arrays used for energy and derivatives. */

  if ( gb==2 || gb==5 ) psi = vector(0, *natom);
  reff = vector(0, *natom);
  sumdeijda = vector(0, *natom);
  for (i = 0; i < *natom; i++) {
    sumdeijda[i] = 0.0;
    sumepol[i] = 0.0;
    sumelec[i] = 0.0;
    sumevdw[i] = 0.0;
  }

  count = 0;
  totsasa = 0.0;

  /* 
   * Get the "effective" Born radii via the approximate pairwise method
   * Use Eqs 9-11 of Hawkins, Cramer, Truhlar, J. Phys. Chem. 100:19824
   * (1996).
   *
   * The computation is not executed in parallel for the case of gbsa==1
   * due to the loop-carried dependence of the 'count' variable.
   */

  if (gb_debug)
    printf("Effective Born radii:\n");

  /* Loop over all atoms i. */

  for (i = 0; i < *natom; i++) {

    xi = x[3 * i];
    yi = x[3 * i + 1];
    zi = x[3 * i + 2];
    ri = rborn[i] - BOFFSET;
    ri1i = 1. / ri;
    sumi = 0.;
    
    icount2 = 0;

    for (k = lprstr[i]; k < lprstr[i] + lpairs[i] + upairs[i]; k++) {
      
      j = pairlist[k];
      
      xij = xi - x[3 * j];
      yij = yi - x[3 * j + 1];
      zij = zi - x[3 * j + 2];
      r2 = xij * xij + yij * yij + zij * zij;
      if ( r2 > rgbmaxpsmax2 ) continue;
      sj = fs[j] * (rborn[j] - BOFFSET);
      if (r2 > (rgbmax+sj)*(rgbmax+sj)) continue;
      nb1[icount2] = j;
      r2x[icount2] = r2;
      vectmp1[icount2] = sj;
      icount2 ++;

      if (gbsa == 1) {
	if ((P0[i] + P0[j])*(P0[i] + P0[j]) > r2) {
	  if (P0[i] > 2.5 && P0[j] > 2.5) {

	    ineighbor[count] = j + 1;
	    /* this +1 is VERY important */
	    count += 1;
	  }
	}
      } /*end gbsa-if */


    }
    for (k=0;k<icount2;k++) vrinv[k] = 1./sqrt(r2x[k]);  /* inverse dists     */
    
    k1=k2=k3=k4=k5=0;
    for (k=0;k<icount2;k++)
      {
	dij = r2x[k]*vrinv[k];
	sj = vectmp1[k];
	if ( dij > rgbmax -sj) 
	  { 
	    ivd1[k1] = vrinv[k];              /* dij1i                        */
	    d1[k1] = dij;                     /* dij                          */
	    vsj1[k1] = sj;                    /*  sj                          */
	    vuij[k1] = dij - vsj1[k1];        /* 1/uij to be inverted later   */
	    vlog1[k1] = vuij[k1];             /* store argument for logarithm */
	    k1++;
	  }
	else if (dij > 4.0* sj)
	  {
	    ivd2[k2] = vrinv[k];              /* dij1i                        */
	    d2[k2] = dij;                     /* dij                          */
	    vsj2[k2] = sj;                    /* sj                           */
	    k2++;
	  }
	else if ( dij > ri + sj)
	  {
	    ivd3[k3] = vrinv[k];
	    vsj3[k3] = sj;
	    k3++;
	  }
	else if (dij > fabs(ri -sj))
	  {
	    ivd4[k4] = vrinv[k];
	    vsj4[k4] = sj;
	    k4++;
	  }
	else if (ri<sj)
	  {
	    ivd5[k5] = vrinv[k];
	    vsj5[k5] = sj;
	    k5++;
	  }
      }

    /* compute the "rgbmax" log() in a vectorizable loops plus the "correct" uij's */
    /* any additional modifications seem to yield no further speedup               */
    /* note that due to the small number of short distance pairs, making the rest  */
    /* of the logarithms etc. vectorizable also yields no further speedup for this */
    /* part of the code                                                            */
    
    for (k=0;k<k1;k++) vlog1[k] = log(vlog1[k]);
    for (k=0;k<k1;k++) vlog1[k] += logrgbmax1i;
    for (k=0;k<k1;k++) vuij[k] = 1./vuij[k];

    for (k=0;k<k1;k++) {
      sj = vsj1[k];
      sj2 = sj*sj;
      dij1i = ivd1[k];
      dij = d1[k];
      r2 = dij*dij;
      uij = vuij[k];
      sumi -= 0.125 * dij1i * (1.0 + 2.0 * (dij *uij + vlog1[k]) +
			       rgbmax2i * (r2 - 4.0 * rgbmax * dij - sj2));
			       }

    
    for (k=0;k<k2;k++) {
      dij2i = ivd2[k]*ivd2[k];
      sj = vsj2[k];
      sj2 = sj*sj;
      tmpsd = sj2*dij2i;
      dumbo = TA+tmpsd*(TB+tmpsd*(TC+tmpsd*(TD+tmpsd*TDD)));
      sumi -= sj*tmpsd*dij2i*dumbo;
    }


    for (k=0;k<k3;k++)
      {
	sj = vsj3[k];
	sj2 = sj*sj;
	dij1i = ivd3[k];
	dij = 1./dij1i;
	r2 = dij*dij;
	sumi -= 0.5 * (sj / (r2 - sj2) +
		       0.5 * dij1i * log((dij - sj) / (dij + sj)));

      }

    for (k=0;k<k4;k++)
      {
	dij1i = ivd4[k];
	dij = 1./dij1i;
	r2 = dij*dij;
	sj = vsj4[k];
	sj2 = sj*sj;
	uij = 1. / (dij + sj);
	theta = 0.5 * ri1i * dij1i * (r2 + ri * ri - sj2);
	sumi -= 0.25 * (ri1i * (2. - theta) - uij +
			dij1i * log(ri * uij));
	
      }

    /* part "5" - the total overlap case - totally useless for currenty radii sets */
    for (k=0; k<k5;k++)
      {
	dij1i = ivd5[k];
	dij = 1./dij1i;
	r2 = dij*dij;
	sj = vsj5[k];
	sj2 = sj*sj;
	sumi -= 0.5 * (sj / (r2 - sj2) + 2. * ri1i +
                       0.5 * dij1i * log((sj - dij) / (sj + dij)));
      }
    
    /* Select atom j from the pair list. */

    if( gb==1 ){

      /* "standard" (HCT) effective radii:  */
      reff[i] = 1.0/(ri1i + sumi);
      if( reff[i] < 0.0 ) reff[i] = 30.0;

    } else {

      /* "gbao" formulas:  */

      psi[i] = -ri*sumi;
      reff[i] = 1.0 / (ri1i - tanh( (gbalpha - gbbeta*psi[i] + 
				     gbgamma*psi[i]*psi[i]) *psi[i])/rborn[i] );
    }

    if (gb_debug)
      printf("%d\t%15.7f\t%15.7f\n", i + 1, rborn[i], reff[i]);

    if (gbsa == 1) {
      ineighbor[count] = 0;
      count = count + 1;
    }
  }

  t2 = second();
  teffr += t2 - t1;
  t1 = t2;

  /*
   * Main LCPO stuff follows.  No parallelization using OpenMP due to
   * the loop-carried dependence of the 'count variable.
   */

  if (gbsa == 1) {
    totsasa = 0.0;
    count = 0;
    for (i = 0; i < *natom; i++) {
      if (ineighbor[count] == 0) {
	count = count + 1;
      } else {

	si = PIx4 * P0[i] * P0[i];
	sumAij = 0.0;
	sumAjk = 0.0;
	sumAjk2 = 0.0;
	sumAijAjk = 0.0;
	sumdAijddijdxi = 0.0;
	sumdAijddijdyi = 0.0;
	sumdAijddijdzi = 0.0;
	sumdAijddijdxiAjk = 0.0;
	sumdAijddijdyiAjk = 0.0;
	sumdAijddijdziAjk = 0.0;

	icount = count;

      L70:j = ineighbor[count] - 1;
				/* mind the -1 , fortran again */
      xi = x[3 * i];
      yi = x[3 * i + 1];
      zi = x[3 * i + 2];

      xj = x[3 * j];
      yj = x[3 * j + 1];
      zj = x[3 * j + 2];



      r2 = (xi - xj) * (xi - xj) + (yi - yj) * (yi - yj) + (zi -
							    zj) *
	(zi - zj);
      dij1i = 1. / sqrt(r2);
      rij = r2 * dij1i;
      tmpaij =
	P0[i] - rij * 0.5 - (P0[i] * P0[i] -
			     P0[j] * P0[j]) * 0.5 * dij1i;

      Aij = PIx2 * P0[i] * tmpaij;

      dAijddij =
	PIx1 * P0[i] * (dij1i * dij1i *
			(P0[i] * P0[i] - P0[j] * P0[j]) - 1.0);

      dAijddijdxj = dAijddij * (xj - xi) * dij1i;
      dAijddijdyj = dAijddij * (yj - yi) * dij1i;
      dAijddijdzj = dAijddij * (zj - zi) * dij1i;

      sumAij = sumAij + Aij;

      count2 = icount;
      sumAjk2 = 0.0;
      sumdAjkddjkdxj = 0.0;
      sumdAjkddjkdyj = 0.0;
      sumdAjkddjkdzj = 0.0;

      p3p4Aij = -surften * (P3[i] + P4[i] * Aij);

      L80:k = ineighbor[count2] - 1;
				/*same as above, -1 ! */
      if (j == k)
	goto L85;

      xk = x[3 * k];
      yk = x[3 * k + 1];
      zk = x[3 * k + 2];

      rjk2 =
	(xj - xk) * (xj - xk) + (yj - yk) * (yj - yk) + (zj -
							 zk) *
	(zj - zk);

      djk1i = 1.0 / sqrt(rjk2);

      rjk = rjk2 * djk1i;

      if (P0[j] + P0[k] > rjk) {
	vdw2dif = P0[j] * P0[j] - P0[k] * P0[k];
	tmpajk = 2.0 * P0[j] - rjk - vdw2dif * djk1i;

	Ajk = PIx1 * P0[j] * tmpajk;

	sumAjk = sumAjk + Ajk;
	sumAjk2 = sumAjk2 + Ajk;

	dAjkddjk =
	  PIx1 * P0[j] * djk1i * (djk1i * djk1i * vdw2dif -
				  1.0);

	dAjkddjkdxj = dAjkddjk * (xj - xk);
	dAjkddjkdyj = dAjkddjk * (yj - yk);
	dAjkddjkdzj = dAjkddjk * (zj - zk);

	f[3 * k] = f[3 * k] + dAjkddjkdxj * p3p4Aij;
	f[3 * k + 1] = f[3 * k + 1] + dAjkddjkdyj * p3p4Aij;
	f[3 * k + 2] = f[3 * k + 2] + dAjkddjkdzj * p3p4Aij;

	sumdAjkddjkdxj = sumdAjkddjkdxj + dAjkddjkdxj;
	sumdAjkddjkdyj = sumdAjkddjkdyj + dAjkddjkdyj;
	sumdAjkddjkdzj = sumdAjkddjkdzj + dAjkddjkdzj;

      }

      L85:count2 = count2 + 1;
      if (ineighbor[count2] != 0) {
	goto L80;
      } else {
	count2 = icount;
      }


      sumAijAjk = sumAijAjk + Aij * sumAjk2;

      sumdAijddijdxi = sumdAijddijdxi - dAijddijdxj;
      sumdAijddijdyi = sumdAijddijdyi - dAijddijdyj;
      sumdAijddijdzi = sumdAijddijdzi - dAijddijdzj;
      sumdAijddijdxiAjk =
	sumdAijddijdxiAjk - dAijddijdxj * sumAjk2;
      sumdAijddijdyiAjk =
	sumdAijddijdyiAjk - dAijddijdyj * sumAjk2;
      sumdAijddijdziAjk =
	sumdAijddijdziAjk - dAijddijdzj * sumAjk2;

      lastxj = dAijddijdxj * sumAjk2 + Aij * sumdAjkddjkdxj;
      lastyj = dAijddijdyj * sumAjk2 + Aij * sumdAjkddjkdyj;
      lastzj = dAijddijdzj * sumAjk2 + Aij * sumdAjkddjkdzj;

      dAidxj =
	surften * (P2[i] * dAijddijdxj +
		   P3[i] * sumdAjkddjkdxj + P4[i] * lastxj);
      dAidyj =
	surften * (P2[i] * dAijddijdyj +
		   P3[i] * sumdAjkddjkdyj + P4[i] * lastyj);
      dAidzj =
	surften * (P2[i] * dAijddijdzj +
		   P3[i] * sumdAjkddjkdzj + P4[i] * lastzj);

      f[3 * j] = f[3 * j] + dAidxj;
      f[3 * j + 1] = f[3 * j + 1] + dAidyj;
      f[3 * j + 2] = f[3 * j + 2] + dAidzj;

      count = count + 1;
      if (ineighbor[count] != 0) {
	goto L70;
      } else {
	count = count + 1;
      }

      Ai = P1[i] * si + P2[i] * sumAij + P3[i] * sumAjk +
	P4[i] * sumAijAjk;

      dAidxi =
	surften * (P2[i] * sumdAijddijdxi +
		   P4[i] * sumdAijddijdxiAjk);
      dAidyi =
	surften * (P2[i] * sumdAijddijdyi +
		   P4[i] * sumdAijddijdyiAjk);
      dAidzi =
	surften * (P2[i] * sumdAijddijdzi +
		   P4[i] * sumdAijddijdziAjk);

      f[3 * i] = f[3 * i] + dAidxi;
      f[3 * i + 1] = f[3 * i + 1] + dAidyi;
      f[3 * i + 2] = f[3 * i + 2] + dAidzi;

      totsasa = totsasa + Ai;
				/*  printf("totsasa up to %d %f\n",i,totsasa); */
      }
    }
    /*      printf("SASA %f , ESURF %f \n",totsasa,totsasa*surften); */
    esurf = totsasa * surften;
  }

  /*
   * Compute the GB, Coulomb and Lennard-Jones energies and derivatives.
   * When in multi-threaded mode in this nest of loops sum data relative only
   * To atom i so that OpenMP can be used to parallelize the outer loop.
   * Perform excluded list management when in single-threaded mode only.
   */

  iexcl = 0;



  {

    /*
     * Allocate and initialize the iexw array used for skipping excluded
     * atoms.  Note that because of the manner in which iexw is used, it
     * is necessary to initialize it before only the first iteration of
     * the following loop.
     */

    iexw = ivector(-1, *natom);
    for (i = -1; i < *natom; i++) {
      iexw[i] = -1;
    }



    /* Loop over all atoms i. */

    for (i = 0; i < *natom; i++) {

      ri = reff[i];
      qi = q[i];


      if (!frozen[i]) {
	expmkf = exp(-KSCALE * (*kappa) * ri) / (*diel_ext);
	dielfac = 1.0 - expmkf;
	qi2h = 0.5 * qi * qi;
	qid2h = qi2h * dielfac;
	sumepol[i] += -qid2h / ri;
	sumdeijda[i] += qid2h - KSCALE * (*kappa) * qi2h * expmkf *ri;
      }

      /*
       * Skip the pair calculations if there are no atoms j on the
       * pair list of atom i.
       */

      nexcl = prm->Iblo[i];
      npairs = upairs[i];
      if ( npairs > 0 ) {

	i3 = 3 * i;

	xi = x[i3    ];
	yi = x[i3 + 1];
	zi = x[i3 + 2];

	iaci = prm->Ntypes * (prm->Iac[i] - 1);

	/*
	 * Expand the excluded atom list into the iexw array by storing i
	 * at array address j.  The starting address is either iexcl that
	 * gets updated with each iteration of the i-loop, or Istr[i] that
	 * has been created in the mme_init function to permit parallelization.
	 */


	jexcl = iexcl;


	jexcl_last = jexcl + nexcl;
	for (j = jexcl; j < jexcl_last; j++) {
	  iexw[prm->ExclAt[j] - 1] = i;
	}

	epol = elec = evdw = 0.0;

	/* Select atoms j from the pair list. */
	icount2 = 0;
	for (k = uprstr[i]; k < uprstr[i] + npairs; k++) {

	  j = pairlist[k];

	  j3 = 3 * j;

	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];
	  r2 = xij * xij + yij * yij + zij * zij;
	  rj = reff[j];
	  rb2 = ri*rj;
	  vfgbi[icount2] = rb2;
	  vefac[icount2] = -0.25*r2/rb2;
	  r2x[icount2] = r2;
	  nb1[icount2] = j;
	  vectmp2[icount2] = qi*q[j];
	  vectmp3[icount2] = ri *rj + 0.25* r2;
	  icount2++;
	}

	/* Do most of the nbond - calculation in a vectorizable way. This is mostly useless except */
	/* for the exp()'s and (inv)sqrt()'s                                                       */
		for (k=0; k < icount2 ; k++) vefac[k] = exp(vefac[k]); /* efac */
		for (k=0; k < icount2 ; k++) vrinv[k] = 1./sqrt(r2x[k]); /* inverse distances */
		for (k=0; k < icount2 ; k++) vfgbi[k] = vfgbi[k]*vefac[k] + r2x[k]; 
		for (k=0; k < icount2 ; k++) vfgbi[k] = 1./sqrt(vfgbi[k]); /* vfgbi complete */
		for (k=0; k < icount2 ; k++) vexpmkf[k] = kkscale/vfgbi[k];
		for (k=0; k < icount2 ; k++) vectmp1[k] = 1.0 - vexpmkf[k];
		for (k=0; k < icount2 ; k++) vexpmkf[k] = exp(vexpmkf[k]) * invdiel_ext; /* expmkf complete */
		for (k=0; k < icount2 ; k++) vectmp1[k] = 1.0 - vexpmkf[k]* vectmp1[k];
		for (k=0; k < icount2 ; k++) vectmp1[k] = vectmp2[k]*vectmp1[k]*vfgbi[k]*vfgbi[k]*vfgbi[k]; /* now contains temp6*/
		for (k=0; k < icount2 ; k++) vectmp7[k] = -vectmp2[k]*(1.0 - vexpmkf[k])*vfgbi[k]; /*delta epol */
		for (k=0; k < icount2 ; k++) vectmp4[k] =  vectmp1[k] * (1.0 - 0.25 * vefac[k]); /* delta de_pol */
		for (k=0; k < icount2 ; k++) vectmp3[k] = 0.5 * vefac[k] * vectmp1[k] *vectmp3[k]; /* vectmp3 gets temp5 */
		
	  
	icount2 = 0;		
	for (k = uprstr[i]; k < uprstr[i] + npairs; k++) {

	  j = pairlist[k];

	  j3 = 3 * j;

	  /* Continue computing the non-diagonal energy term. */

	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];
	  rj = reff[j];
	  qiqj = vectmp2[icount2];

	  epol += vectmp7[icount2];

	  de = vectmp4[icount2];
	  
	  temp5 = vectmp3[icount2];

	  sumdeijda[i] += ri * temp5;


	  sumdeijda[j] += rj * temp5;


	  if (iexw[j] != i) {

	    rinv = vrinv[icount2];
	    r2inv = rinv * rinv;

	    /*  gas-phase Coulomb energy:  */

	    eel = qiqj * rinv;
	    elec += eel;
	    de -= eel * r2inv;

	    /* Lennard-Jones energy:   */

	    ic = prm->Cno[iaci + prm->Iac[j] - 1] - 1;
	    if( ic >= 0 ){
	      r6inv = r2inv * r2inv * r2inv;
	      f6 = prm->Cn2[ic] * r6inv;
	      f12 = prm->Cn1[ic] * r6inv * r6inv;
	      evdw += f12 - f6;
	      de -= (12. * f12 - 6. * f6) * r2inv;
	    }
	  }
	  icount2 ++;
	  /*
	   * When in single-threaded mode sum to the gradient vector the
	   * derivatives of Dij that are computed relative to the cartesian
	   * coordinates of atoms i and j.  When in multi-threaded mode sum
	   * to the gradient vector the derivatives that are computed relative
	   * to the cartesian coordinates of atom i.
	   */

	  dedx = de * xij;
	  dedy = de * yij;
	  dedz = de * zij;

	  f[i3    ] += dedx;
	  f[i3 + 1] += dedy;
	  f[i3 + 2] += dedz;


	  f[j3    ] -= dedx;
	  f[j3 + 1] -= dedy;
	  f[j3 + 2] -= dedz;

	}

	/* Sum the energy terms. */

	sumepol[i] += epol;
	sumelec[i] += elec;
	sumevdw[i] += evdw;
      }

      /*
       * When in single-threaded mode update iexcl.  This update is unnecessary
       * in multi-threaded mode because jexcl is loaded from Istr[i] not iexcl.
       */


      iexcl += nexcl;

    }

    /* Free the iexw array within this potentially parallel region of code. */

    free_ivector(iexw, -1, *natom);

    t2 = second();
    tnonb += t2 - t1;
    t1 = t2;

    /*
     * Compute the derivatives of the effective radius Ri of atom i
     * with respect to the cartesian coordinates of each atom j.
     *
     * When in single-threaded mode sum all of these derivatives into the
     * gradient vector.  When in multi-threaded mode sum into the gradient
     * vector only the derivatives with respect to the cartesian coordinates
     * of atom i (which are in fact the sums of all other derivatives).
     */


    /* Loop over all atoms i. */

    for (i = 0; i < *natom; i++) {

      /*
       * Don't calculate derivatives of the effective radius of atom i
       * if atom i is frozen or if there are no pair atoms j associated
       * with atom i.
       */

      npairs = lpairs[i] + upairs[i];
      if ( frozen[i] || !npairs ) continue;

      i3 = 3 * i;

      xi = x[i3    ];
      yi = x[i3 + 1];
      zi = x[i3 + 2];
      ri = rborn[i] - BOFFSET;
      ri1i = 1. / ri;

      sumda = sumdeijda[i];

      if ( gb>1 ) {

	ri = rborn[i] - BOFFSET;
	thi = tanh( (gbalpha - gbbeta*psi[i] + gbgamma*psi[i]*psi[i])*psi[i] );
	sumda *= (gbalpha -2.0*gbbeta*psi[i] + 3.0*gbgamma*psi[i]*psi[i])
	  *(1.0 - thi*thi)*ri/rborn[i];
      }

      /* Initialize the derivative accumulators. */

      daix = daiy = daiz = 0.0;

      /* Select atom j from the pair list. */
      icount2=0;
      for (k = lprstr[i]; k < lprstr[i] + npairs; k++) {

	j = pairlist[k];

	j3 = 3 * j;

	xij = xi - x[j3    ];
	yij = yi - x[j3 + 1];
	zij = zi - x[j3 + 2];
	r2 = xij * xij + yij * yij + zij * zij;
	if ( r2 > rgbmaxpsmax2 ) continue;
	sj = fs[j] * (rborn[j] - BOFFSET);
	if (r2 > (rgbmax+sj)*(rgbmax+sj)) continue;
	r2x[icount2] = r2;
	nb1[icount2] = j; /* all neighbor indices of atom i stored here */
	vsj[icount2] = sj;
	icount2++;
      }
      
      for (k=0; k<icount2;k++) vrinv[k] = 1./sqrt(r2x[k]); /* the inverse distances again  */

      /* what follows is analogous to the radii calculation  */
      /* principle: first generate lists of what is to be done, then loop over those lists */
      
      k1 = k2 = k3 = k4 = k5 = 0;
      for (k=0;k<icount2;k++)
	{
	  j=nb1[k];
	  r2=r2x[k];
	  sj = vsj[k];
	  sj2 = sj*sj;
	  dij1i = vrinv[k];
	  dij = r2*dij1i;

	  if (dij > rgbmax -sj)
	    {
	      nbl1[k1] = j;
	      ivd1[k1] = dij1i;
	      vsj1[k1] = sj;
	      vectmp1[k1] = dij - sj;
	      vectmp3[k1] = r2;
	      k1++;
	    }
	  else if (dij > 4.0*sj)
	    {
	      nbl2[k2] = j;
	      vsj2[k2] = sj;
	      ivd2[k2] = dij1i*dij1i;
	      vectmp5[k2] = sj2*dij1i*dij1i;
	      k2++;
	    }
	  else if (dij >ri + sj)
	    {
	      nbl3[k3] = j;
	      vsj3[k3] = sj;
	      ivd3[k3] = vrinv[k];
	      k3++;
	    }
	  else if (dij > fabs(ri -sj))
	    {
	      nbl4[k4] = j;
	      k4++;
	    }
	  else if ( ri <sj)
	    {
	      nbl5[k5] = j;
	      k5++;
	    }
	}

      /* after the following loop, vectmp4 contains the "datmps" for all  */
      /* neighbors j that lie within the rgbmax-smoothing region        */
      for (k=0;k<k1;k++) vectmp1[k] = 1./vectmp1[k];
      for (k=0;k<k1;k++) vectmp2[k] = log(vectmp1[k]);
      for (k=0;k<k1;k++) vectmp2[k] -= logrgbmax1i;
      for (k=0;k<k1;k++) vectmp2[k] *= 2.0;
      for (k=0;k<k1;k++) vectmp4[k] = vectmp1[k]*vectmp1[k] - rgbmax2i;
      for (k=0;k<k1;k++) vectmp4[k] = vectmp4[k]*(vectmp3[k] + vsj1[k]*vsj1[k]);
      for (k=0;k<k1;k++) vectmp4[k] -= vectmp2[k];
      for (k=0;k<k1;k++) vectmp4[k] = 0.125*vectmp4[k]*ivd1[k]*ivd1[k]*ivd1[k];

      for (k=0;k<k1;k++)
	{
	  j=nbl1[k];
	  j3=3*j;
	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];

	  datmp = vectmp4[k];
	  
	  daix += xij * datmp;
	  daiy += yij * datmp;
	  daiz += zij * datmp;

	  datmp *= sumda;
	  f[j3    ] += xij * datmp;
	  f[j3 + 1] += yij * datmp;
	  f[j3 + 2] += zij * datmp;
	}

      /* the next three lines calculate datmp for all neightbors within the                 */
      /*"Taylor approximated"   region. vectmp7 contains the these datmps                   */
      for (k=0;k<k2;k++) vectmp6[k] =  TE + vectmp5[k]*(TF+vectmp5[k]*(TG+vectmp5[k]*(TH+vectmp5[k]*THH)));
      for (k=0;k<k2;k++) vectmp7[k] = vectmp5[k]*vectmp6[k];
      for (k=0;k<k2;k++) vectmp7[k] = vectmp7[k]*vsj2[k]*ivd2[k]*ivd2[k];
      

      
      for (k=0;k<k2;k++)
	{
	  j=nbl2[k];
	  j3=3*j;
	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];

	  datmp = vectmp7[k];
	  
	  daix += xij * datmp;
	  daiy += yij * datmp;
	  daiz += zij * datmp;

	  datmp *= sumda;

	  f[j3    ] += xij * datmp;
	  f[j3 + 1] += yij * datmp;
	  f[j3 + 2] += zij * datmp;
	}

      /* for the rest, additional vectorizability-attempts again don't help */
      for (k=0;k<k3;k++)
	{
	  j=nbl3[k];
	  j3=3*j;
	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];
	  r2 = xij*xij + yij*yij + zij*zij;
	  dij1i = ivd3[k];
	  dij = r2*dij1i;
	  dij2i = dij1i*dij1i;
	  dij3i = dij1i*dij2i;
	  sj = vsj3[k];
	  sj2 = sj*sj;

	  temp1 = 1. / (r2 - sj2);
	  datmp = temp1 * sj * (-0.5 * dij2i + temp1)
	    + 0.25 * dij3i* log((dij - sj) / (dij + sj));


	  daix += xij * datmp;
	  daiy += yij * datmp;
	  daiz += zij * datmp;

	  datmp *= sumda;
	  f[j3    ] += xij * datmp;
	  f[j3 + 1] += yij * datmp;
	  f[j3 + 2] += zij * datmp;
	}

      for (k=0;k<k4;k++)
	{
	  j=nbl4[k];
	  j3=3*j;
	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];
	  r2 = xij*xij + yij*yij + zij*zij;
	  dij1i = 1./sqrt(r2);
	  dij = r2*dij1i;
	  dij2i = dij1i*dij1i;
	  dij3i = dij1i*dij2i;
	  sj = fs[j] * (rborn[j] - BOFFSET);
	  sj2 = sj*sj;

	  
	  temp1 = 1. / (dij + sj);

	  datmp = -0.25 * (-0.5 * (r2 - ri * ri + sj2) * dij3i * ri1i *
			 ri1i + dij1i * temp1 * (temp1 - dij1i)
			   - dij3i * log(ri * temp1));

	  daix += xij * datmp;
	  daiy += yij * datmp;
	  daiz += zij * datmp;

	  datmp *= sumda;
	  f[j3    ] += xij * datmp;
	  f[j3 + 1] += yij * datmp;
	  f[j3 + 2] += zij * datmp;
	}
      for (k=0;k<k5;k++)
	{
	  j=nbl5[k];
	  j3=3*j;
	  xij = xi - x[j3    ];
	  yij = yi - x[j3 + 1];
	  zij = zi - x[j3 + 2];
	  r2 = xij*xij + yij*yij + zij*zij;
	  dij1i = 1./sqrt(r2);
	  dij = r2*dij1i;
	  dij2i = dij1i*dij1i;
	  dij3i = dij1i*dij2i;
	  sj = fs[j] * (rborn[j] - BOFFSET);
	  sj2 = sj*sj;

	  	  
	  temp1 = 1. / (dij + sj);
	  
	  datmp = -0.25 * (-0.5 * (r2 - ri * ri + sj2) * dij3i * ri1i *
			   ri1i + dij1i * temp1 * (temp1 - dij1i)
			   - dij3i * log(ri * temp1));
	  
	  daix += xij * datmp;
	  daiy += yij * datmp;
	  daiz += zij * datmp;

	  datmp *= sumda;
	  f[j3    ] += xij * datmp;
	  f[j3 + 1] += yij * datmp;
	  f[j3 + 2] += zij * datmp;
	}
      f[i3    ] -= sumda * daix;
      f[i3 + 1] -= sumda * daiy;
      f[i3 + 2] -= sumda * daiz;
    }

  }

  /* Perform non-logarithmic reduction of the energy terms. */

  epol = elec = evdw = 0.0;

  for (i = 0; i < *natom; i++) {
    epol += sumepol[i];
    elec += sumelec[i];
    evdw += sumevdw[i];
  }

  /* Return elec and evdw through the reference parameters eelt and enb. */

  *eelt = elec;
  *enb = evdw;

  if (gb_debug > 1) {
    printf("gb    epol= %10.3f\n", epol);
    for (k = 0; k < *natom; k++) {
      printf("%15.7f\n", sumdeijda[k]);
    }
    for (k = 0; k < *natom * 3; k += 3) {
      printf("%15.7f%15.7f%15.7f\n", f[k], f[k + 1], f[k + 2]);
    }
  }

  free_vector(reff, 0, *natom);
  free_vector(sumepol, 0, *natom);
  free_vector(sumelec, 0, *natom);
  free_vector(sumevdw, 0, *natom);
  free_vector(sumdeijda, 0, *natom);
  if( gb == 2 || gb==5 ) free_vector(psi, 0, *natom);

  free_vector (vfgbi, 0, *natom);
  free_vector (vefac, 0,*natom);
  free_vector (r2x, 0, *natom);
  free_vector  (vrinv, 0, *natom);
  free_vector (vexpmkf, 0, *natom);
  free_vector  (vectmp1, 0, *natom);
  free_vector  (vectmp2, 0, *natom);
  free_vector  (vectmp3, 0, *natom);
  free_vector  (vectmp4, 0, *natom);
  free_vector  (vectmp5, 0, *natom);
  free_vector  (vectmp6, 0, *natom);
  free_vector  (vectmp7, 0, *natom);
  
  free_ivector (nb1 , 0 , *natom);
  free_ivector (nbl1 , 0 , *natom);
  free_ivector (nbl2 , 0 , *natom);
  free_ivector (nbl3 , 0 , 150);
  free_ivector (nbl4 , 0 , 50);
  free_ivector (nbl5 , 0 , 5);

 /*  
  REAL_T *ivd1,*ivd2,*ivd3,*ivd4,*ivd5;
  REAL_T *vsj1,*vsj2,*vsj3,*vsj4,*vsj5;
  REAL_T *d1,*d2,*vuij,*vlog1;
 */
  free_vector(ivd1,0, *natom);
  free_vector(ivd2,0, *natom);
  free_vector(ivd3,0, 150);
  free_vector(ivd4,0, 50);
  free_vector(ivd5,0, 5);

  free_vector(vsj,0, *natom);
  free_vector(vsj1,0, *natom);
  free_vector(vsj2,0, *natom);
  free_vector(vsj3,0, 150);
  free_vector(vsj4,0, 50);
  free_vector(vsj5,0, 5);
  
  free_vector(d1,0,*natom);
  free_vector(d2,0, *natom);
  free_vector (vuij, 0, *natom);
  free_vector(vlog1, 0, *natom);


  t2 = second();
  tdiag += t2 - t1;
  t1 = t2;

  return (epol);
}