double AmberConformation::solventEnergy(){
  double solvEnergy = 0;

  //constants... may put in .h file
  double offset = 0.09;
  double rgbmax = 25;
  double rgbmax2 = rgbmax*rgbmax;
  double rgbmax1i = 1/rgbmax;
  double rgbmax2i = rgbmax1i*rgbmax1i;
  double rgbmaxpsmax2;

  double cut =  90000.0;
  double cut_inner = 64.0;  

  double extdiel = 78.5;
  double intdiel = 1.0;
  double extdieli = 1.0/extdiel;
  double intdieli = 1.0/intdiel;

  double saltcon = 0.2;
  double kappa = sqrt( 0.10806 * saltcon );
  double surften = 0.005;
  kappa = 0.73*kappa;


  
  //does this need to be so big???
  //int * ineighbor = new int[anumAtoms*anumAtoms];

  for(int i = 0; i < anumAtoms; i++){
    reff[i] = 0;
    sumdeijda[i] = 0;
    //ineighbor[i] = 0;
    //fs[i] = 0;
    //jj[i] = 0;
    // r2x[i] = 0;
//     rjx[i] = 0;
//     vectmp1[i] = 0;
//     vectmp2[i] = 0;
//     vectmp3[i] = 0;
//     vectmp4[i] = 0;
//     vectmp5[i] = 0;
//     psi[i] = 0;
//     j_vals[i] = 0;
//     k_vals[i] = 0;
  }

  //double * fs = new double[anumAtoms];
  double fsmax = 0.0;

  for(int i = 0; i < anumAtoms; i++){
    fs[i] = atoms[i].screen * (atoms[i].radii - offset);
    if(fs[i] > fsmax)
      fsmax = fs[i];
  }

  rgbmaxpsmax2 = (rgbmax+fsmax);
  rgbmaxpsmax2 *= rgbmaxpsmax2;
    
  double xi, yi, zi, xij, yij, zij, ri, rj, rij, ri1i, rj1i, si, si2, reff_i, r2, dij1i, dij2i, dij3i, dij, sj, sj2, sij, uij;
  double de, fgbi, rinv, qi, qiqj, r2inv, rb2, fgbk, expmkf, e, dl, qi2h, qid2h;
  double temp1, temp4, temp5, temp6, temp7, datmp;
  double thi, thi2;
  double tmpsd, dumbo, theta;

  //LCPO variables
  double sumaij, sumajk, sumajk2, sumaijajk, sumdx, sumdy, sumdz, sumidx, sumidy, sumidz, sumkdx, sumkdy, sumkdz, dajkddjk;
  double daijddij, ddx, ddy, ddz, dajx, dajy, dajz, daidxj, daidyj, daidzj, daidxi, daidyi, daidzi, ai;
  double lastxj, lastyj, lastzj;
  
  int icount, kk1, kk2;

  int count = 0;
 
  double ta = 1.0/3.0;
  double tb = 2.0/5.0;
  double tc = 3.0/7.0;
  double td = 4.0/9.0;
  double tdd = 5.0/11.0;
  double te = 4.0/3.0;
  double tf = 12.0/5.0;
  double tg = (3.0*8.0)/7.0;
  double th = (5.0*8.0)/9.0;
  double thh = (3.0*4.0*5.0/11.0);

  double gbalpha = 1.0;
  double gbbeta = 0.8;
  double gbgamma = 4.851;

  double epol = 0.0;
  double elet = 0.0;
  double evdw = 0.0;
  double esurf = 0.0;
  double totsasa = 0.0;

  for(int i = 0; i < anumAtoms; i++){
    //printf("i: %d\n", i);
    xi = coords[i].x;
    yi = coords[i].y;
    zi = coords[i].z;

    reff_i = reff[i];
    ri = atoms[i].radii - offset;
    ri1i = 1/ri;

    si = fs[i];
    si2 = si*si;

    icount = 0;

    for(int j = i+1; j < anumAtoms; j++){
      xij = xi - coords[j].x;
      yij = yi - coords[j].y;
      zij = zi - coords[j].z;

      r2 = xij*xij + yij*yij + zij*zij;
      if(r2 > rgbmaxpsmax2)
	continue;

      jj[icount] = j;
      r2x[icount] = r2;

      icount++;
    }

    kk1 = 0;
    kk2 = 0;

    for(int k = 0; k < icount; k++){
      vectmp1[k] = 1/sqrt(r2x[k]);

      int j = jj[k];
      r2 = r2x[k];
      sj = fs[j];

      dij1i = vectmp1[k];
      dij = r2*dij1i;

      if(dij > rgbmax + si && dij > rgbmax+sj)
	continue;

      rj = atoms[j].radii - offset;

	
      if(dij <= 4*sj){
	vectmp2[kk1] = dij + sj;
	if(dij > ri+sj)
	  vectmp4[kk1] = dij - sj;
	else if(dij > fabs(ri-sj))
	  vectmp4[kk1] = ri;
	else if(ri < sj)
	  vectmp4[kk1] = sj - dij;
	else
	  vectmp4[kk1] = 1;


	kk1++;
      }

      if(dij <= 4*si){
	vectmp3[kk2] = dij + si;
	if(dij > rj+si)
	  vectmp5[kk2] = dij - si;
	else if(dij > fabs(rj-si))
	  vectmp5[kk2] = rj;
	else if(rj < si)
	  vectmp5[kk2] = si - dij;
	else
	  vectmp5[kk2] = 1;


	kk2++;
      }
    }


    for(int k = 0; k < kk1; k++)
      vectmp2[k] = 1/vectmp2[k];

    for(int k = 0; k < kk2; k++)
      vectmp3[k] = 1/vectmp3[k];

    for(int k = 0; k < kk1; k++){
      vectmp4[k] = vectmp2[k]*vectmp4[k];
      vectmp4[k] = log(vectmp4[k]);
    }

    for(int k = 0; k < kk2; k++){
      vectmp5[k] = vectmp3[k]*vectmp5[k];
      vectmp5[k] = log(vectmp5[k]);
    }

    kk1 = 0;
    kk2 = 0;

    for(int k = 0; k < icount; k++){
      int j = jj[k];
      r2 = r2x[k];

      rj = atoms[j].radii - offset;
      rj1i = 1/rj;
      sj = fs[j];

      sj2 = sj * sj;

      xij = xi - coords[j].x;
      yij = yi - coords[j].y;
      zij = zi - coords[j].z;

      dij1i = vectmp1[k];
      dij = r2*dij1i;

      if(dij <= rgbmax + sj){
	//printf("dij: %8.10f\n", dij);
	if(dij > (rgbmax - sj)){
	  uij = 1.0/(dij - sj);
	  //printf("uij: %8.10f\n", uij);
	  reff_i = reff_i - 0.125 * dij1i * (1 + 2 * dij * uij + 
					     rgbmax2i * (r2 - 4 * rgbmax * dij - sj2) +
					     2 * log((dij-sj)*rgbmax1i));

	  //printf("reff_i: %8.15f\n", reff_i);
	}

	else if(dij > 4*sj){
	  dij2i = dij1i*dij1i;
	  tmpsd = sj2*dij2i;
	  dumbo = ta+tmpsd*(tb+tmpsd*(tc+tmpsd*(td+tmpsd*tdd)));
	  
	  //printf("dumbo: %8.10f\n", dumbo);
	  reff_i = reff_i - tmpsd*sj*dij2i*dumbo;
	}

	else if(dij > ri+sj){
	  //printf("blah\n");
	  reff_i = reff_i - (1.0/2.0)*(sj/(r2-sj2) + (1.0/2.0)*dij1i*vectmp4[kk1]);
	  kk1++;
	}

	else if(dij > fabs(ri-sj)){
	  //printf("ri1i: %8.10f dij1i: %8.10f r2: %8.10f ri: %8.10f, sj2: %8.10f\n", ri1i, dij1i, r2, ri, sj2);
	  theta = (1.0/2.0)*ri1i*dij1i*(r2 + ri*ri - sj2);
	  //printf("theta: %8.10f\n", theta);
	  reff_i = reff_i - (1.0/4.0)*(ri1i*(2-theta) - vectmp2[kk1] + dij1i*vectmp4[kk1]);
	  kk1++;
	}
	
	else if(ri < sj){
	  //printf("blah2\n");
	  reff_i = reff_i - (1.0/2.0)*(sj/(r2-sj2) + 2.0*ri1i + (1.0/2.0)*dij1i*vectmp4[kk1]);
	  kk1++;
	}
	else{
	  //printf("blah3\n");
	  kk1++;
	}
	  
	//printf("reff_i: %8.10f\n", reff_i);
      }

      if(dij > rgbmax + si)
	continue;

      if(dij > rgbmax - si){
	uij = 1.0/(dij - si);
	reff[j] = reff[j] - 0.125*dij1i*(1.0 + 2.0*dij*uij +
					 rgbmax2i * (r2 - 4.0*rgbmax*dij - si2) +
					 2.0 * log((dij-si)*rgbmax1i));
      }

      else if(dij > 4.0*si){
	dij2i = dij1i*dij1i;
	tmpsd = si2*dij2i;
	dumbo = ta+tmpsd*(tb+tmpsd*(tc+tmpsd*(td+tmpsd*tdd)));
	reff[j] = reff[j] - tmpsd*si*dij2i*dumbo;
      }

      else if(dij > rj+si){
	reff[j] = reff[j] - (1.0/2.0)*(si/(r2-si2) +
				       (1.0/2.0)*dij1i*vectmp5[kk2]);
	kk2++;
      }
	
      else if(dij > fabs(rj-si)){
	theta = (1.0/2.0)*rj1i*dij1i*(r2 + rj*rj - si2);
	reff[j] = reff[j] - (1.0/4.0)*(rj1i*(2.0-theta) - vectmp3[kk2] + dij1i*vectmp5[kk2]);
	kk2++;
      }

      else if(rj < si){
	reff[j] = reff[j] - (1.0/2.0)*(si/(r2-si2) + 2.0*rj1i
				       + (1.0/2.0)*dij1i*vectmp5[kk2]);
	kk2++;
      }
      
      else
	kk2++;
    
    }
    
    reff[i] = reff_i;
      
  }


  //now do the Onufriev correction:

  for(int i = 0; i < anumAtoms; i++){
    //printf("reff[%i]: %8.15f\n", i, reff[i]);
    ri = atoms[i].radii - offset;
    ri1i = 1.0/ri;
    psi[i] = -ri*reff[i];
    //printf("psii: %8.15f\n", psii);
    //printf("gbalpha: %8.15f gbgamma: %8.15f gbbeta: %8.15f\n", gbalpha, gbgamma, gbbeta);
    //printf("before tanh: %8.15f\n", (gbalpha + gbgamma*psii*psii - gbbeta*psii)*psii);
    reff[i] = ri1i - tanh( (gbalpha + gbgamma*psi[i]*psi[i] - gbbeta*psi[i])*psi[i])/atoms[i].radii;
    reff[i] = 1.0/reff[i];

    //printf("reff[%i]: %8.15f\n", i, reff[i]);
  }
 

  //now calculate the energy


  for(int i = 0; i < anumAtoms; i++){

    

    xi = coords[i].x;
    yi = coords[i].y;
    zi = coords[i].z;
    qi = atoms[i].charge;
    ri = reff[i];
    
    icount = 0;

    for(int j = i+1; j < anumAtoms; j++){
      xij = xi - coords[j].x;
      yij = yi - coords[j].y;
      zij = zi - coords[j].z;
      r2 = xij*xij + yij*yij + zij*zij;
      
      if(r2 > cut)
	continue;

      jj[icount] = j;
      r2x[icount] = r2;
      rjx[icount] = reff[j];
      icount++;
    }

    //printf("icount: %i\n", icount);

    for(int k = 0; k < icount; k++){
      vectmp1[k] = exp(-r2x[k]/(4.0*ri*rjx[k]));
      vectmp3[k] = r2x[k] + rjx[k]*ri*vectmp1[k]; //non sqrt'd fgb
      vectmp2[k] = 1.0/sqrt(vectmp3[k]);
      vectmp3[k] = 1.0/sqrt(r2x[k]);

      if(kappa != 0){
	vectmp4[k] = exp(-kappa/vectmp2[k]);
      }
    }
      

    for(int k = 0; k < icount; k++){
      int j = jj[k];
      de = 0.0;
      r2 = r2x[k];
      rj = rjx[k];
      temp1 = vectmp1[k];
      //printf("temp1: %8.10f\n", temp1);
      fgbi = vectmp2[k];
      //printf("fgbi: %8.10f\n", fgbi);
      rinv = vectmp3[k];
      //printf("rinv: %8.10f\n", rinv);

      xij = xi - coords[j].x;
      yij = yi - coords[j].y;
      zij = zi - coords[j].z;

      //going to ignore frespa, should equal 1 during minimization

      r2inv = rinv*rinv;

      qiqj = qi* atoms[j].charge;
      rb2 = ri*rj;


      if(kappa == 0){
	fgbk = 0;
	expmkf = extdiel;
      }
      else{
	fgbk = -kappa/fgbi;
	expmkf = vectmp4[k]*extdieli;
      }

      dl = intdieli - expmkf;

      e = -qiqj*dl*fgbi;
      //printf("e: %8.10f\n", e);
      epol = epol + e;
      //printf("epol: %8.10f\n", epol);


      temp4 = fgbi*fgbi*fgbi;
      temp6 = -qiqj*temp4*(dl + fgbk*expmkf);
      de = temp6*(1.0 - 0.25*temp1);
      //printf("i: %d j: %d de: %8.15f zij: %8.10f\n", i, j, de, zij);
      forces[j].x -= de*xij;
      forces[j].y -= de*yij;
      forces[j].z -= de*zij;

      forces[i].x += de*xij;
      forces[i].y += de*yij;
      forces[i].z += de*zij;

      //printf("i: %d j: %d de: %8.10f\n", i, j, de);
      temp5 = 0.5*temp1*temp6*(rb2 + 0.25*r2);

      sumdeijda[i] = sumdeijda[i] + ri*temp5;
      //printf("sumdeijda[%i]: %8.15f\n", i, sumdeijda[i]);
      sumdeijda[j] = sumdeijda[j] + rj*temp5;

    }
  
  }
  
  double f_x, f_y, f_z, f_xi, f_yi, f_zi;

  

  for(int i = 0; i < anumAtoms; i++){
    f_xi = f_yi = f_zi = 0.0;
    qi = atoms[i].charge;
    expmkf = exp( -kappa * reff[i])*extdieli;
    dl = intdieli - expmkf;
    qi2h = 0.5*qi*qi;
    qid2h = qi2h*dl;
    epol = epol - qid2h/reff[i];

    temp7 = -sumdeijda[i] + qid2h - kappa*qi2h*expmkf*reff[i];

    xi = coords[i].x;
    yi = coords[i].y;
    zi = coords[i].z;
    ri = atoms[i].radii - offset;
    ri1i = 1/ri;
    
    //printf("psi[%i]: %8.15f\n", i, psi[i]);
    thi = tanh( (gbalpha + gbgamma*psi[i]*psi[i] - gbbeta*psi[i])*psi[i]);
    thi2 = (gbalpha + 3.0*gbgamma*psi[i]*psi[i] - 2.0*gbbeta*psi[i])*(1.0 - thi*thi)*ri/atoms[i].radii;
    //printf("thi2: %8.15f\n", thi2);

    icount = 0;
    for(int j = 0; j < anumAtoms; j++){
      if(i == j)
	continue;

      xij = xi - coords[j].x;
      yij = yi - coords[j].y;
      zij = zi - coords[j].z;

      r2 = xij*xij + yij*yij + zij*zij;
      if(r2 > rgbmaxpsmax2)
	continue;

      jj[icount] = j;
      r2x[icount] = r2;
      vectmp1[icount] = 1.0/sqrt(r2);
      icount++;
    }

    kk1 = 0;

    //printf("icount: %i\n", icount);
    for(int k = 0; k < icount; k++){
      int j = jj[k];
      //printf("1 k: %i j: %i\n", k, jj[k]);
      r2  = r2x[k];

      sj = fs[j];

      dij1i = vectmp1[k];
      dij = r2*dij1i;
      sj2 = sj*sj;

      if(dij > 4.0*sj)
	continue;

      vectmp3[kk1] = dij+sj;

      if(dij > ri+sj){
	vectmp2[kk1] = r2 - sj2;
	vectmp4[kk1] = dij - sj;
      }
      else if (dij > fabs(ri-sj)){
	vectmp2[kk1] = dij + sj;
	vectmp4[kk1] = ri;
      }
      else if(ri < sj){
	vectmp2[kk1] = r2 - sj2;
	vectmp4[kk1] = sj - dij;
      }
      else {
	vectmp2[kk1] = 1.0;
	vectmp4[kk1] = 1.0;
      }
      
      kk1++;
    }

    //printf("kk1: %i\n", kk1);

    for(int k =0; k < kk1; k++){
      vectmp2[k] = 1.0/vectmp2[k];
      vectmp3[k] = 1.0/vectmp3[k];
      vectmp4[k] = log(vectmp4[k] * vectmp3[k]);
    }


    kk1 = 0;
    
    for(int k = 0; k < icount; k++){
      int j = jj[k];
      //printf("i: %d k: %d j: %d\n", i, k, j);
      r2 = r2x[k];
      
      //printf("3 k: %i j: %i\n", k, j);
      xij = xi - coords[j].x;
      yij = yi - coords[j].y;
      zij = zi - coords[j].z;

      dij1i = vectmp1[k];
      //printf("dij1i: %8.15f\n", dij1i);
      dij = r2*dij1i;
      sj = fs[j];
      if(dij > rgbmax+sj)
	continue;

      sj2 = sj*sj;

      dij2i = dij1i*dij1i;
      dij3i = dij2i*dij1i;

      if(dij > rgbmax - sj){
	temp1 = 1.0/(dij-sj);
	datmp = 0.125*dij3i * ((r2 + sj2) * (temp1*temp1 - rgbmax2i) - 2.0*log(rgbmax*temp1));
      }
      else if(dij > 4.0*sj){
	tmpsd = sj2*dij2i;
	dumbo = te+tmpsd*(tf+tmpsd*(tg+tmpsd*(th+tmpsd*thh)));
	datmp = tmpsd*sj*dij2i*dij2i*dumbo;
      }
      else if(dij > ri+sj){
	datmp = vectmp2[kk1]*sj*(-0.5*dij2i + vectmp2[kk1]) + 0.25*dij3i*vectmp4[kk1];
	kk1++;
      }
      else if(dij > fabs(ri-sj)){
	datmp = -0.25*(-0.5*(r2 - ri*ri + sj2)*dij3i*ri1i*ri1i + dij1i*vectmp2[kk1]*(vectmp2[kk1] - dij1i) - dij3i*vectmp4[kk1]);
	kk1++;
      }
      else if(ri < sj){
	datmp = -0.5*(sj*dij2i*vectmp2[kk1] - 2.0*sj*vectmp2[kk1]*vectmp2[kk1] - 0.5*dij3i*vectmp4[kk1]);
	kk1++;
      }
      else{
	datmp = 0.0;
	kk1++;
      }

      //printf("datmp: %8.15f\n", datmp);

      datmp = -datmp*temp7;
      datmp = datmp*thi2;

      f_x = xij*datmp;
      f_y = yij*datmp;
      f_z = zij*datmp;

      //printf("f_x: %8.15f\n", f_x);
      //printf("f_y: %8.15f\n", f_y);
      //printf("f_z: %8.15f\n", f_z);

      forces[j].x += f_x;
      forces[j].y += f_y;
      forces[j].z += f_z;

      f_xi -= f_x;
      f_yi -= f_y;
      f_zi -= f_z;

      //printf("f_xi: %8.15f\n", f_xi);
      //printf("f_yi: %8.15f\n", f_yi);
      //printf("f_zi: %8.15f\n", f_zi);
      
      
      //define the neighbor list for LCPO
      dij = sqrt(r2x[k]);
      if((vdwRad[i] + vdwRad[j]) > dij){
	if(vdwRad[i] > 2.5 && vdwRad[j] > 2.5){
	  ineighbor[count] = j;
	  count++;
	}
      }

    }

    //printf("f_xi: %8.15f\n", f_xi);
    //printf("f_yi: %8.15f\n", f_yi);
    //printf("f_zi: %8.15f\n", f_zi);
    
    forces[i].x += f_xi;
    forces[i].y += f_yi;
    forces[i].z += f_zi;
      
  
    ineighbor[count] = -1;
    count++;
    
    //printf("count: %d\n", count);
  }

  //printf("count: %d\n", count);
  //LCPO method to calculate surface area
  int max_count = count;
  count = 0;
  int num_j_vals, num_k_vals, count2;

  //int * j_vals = new int[anumAtoms];
  //int * k_vals = new int[anumAtoms];
  double tmpaij, aij, ppa;
  double xjk, yjk, zjk, rjk2, djk1i, rjk, tmpajk, ajk, vdw2dif;
  for(int i = 0; i < anumAtoms; i++){
    //printf("ineighbor[%d]: %d\n", count+1, ineighbor[count]+1);
    if(ineighbor[count] == -1){
      count++;
    }
    else{
      si = 4.0*M_PI*vdwRad[i]*vdwRad[i];
      sumaij = 0.0;
      sumajk = 0.0;
      sumajk2 = 0.0;
      sumaijajk = 0.0;
      sumidx = 0.0;
      sumidy = 0.0;
      sumidz = 0.0;
      sumdx = 0.0;
      sumdy = 0.0;
      sumdz = 0.0;

      icount = count;
      num_j_vals = 0;
      while(ineighbor[count] != -1){
	int j = ineighbor[count];
	j_vals[num_j_vals] = j;
	//printf("count: %d j: %d\n", count+1, j+1);
	count++;
	num_j_vals++;
      }
      count++;
      //printf("num_j_vals: %d\n", num_j_vals);

      xi = coords[i].x;
      yi = coords[i].y;
      zi = coords[i].z;

      for(int jjj = 0; jjj < num_j_vals; jjj++){
	int j = j_vals[jjj];
	xij = xi - coords[j].x;
	yij = yi - coords[j].y;
	zij = zi - coords[j].z;

	r2 = xij*xij + yij*yij + zij*zij;
	dij1i = 1.0/sqrt(r2);
	//printf("dij1i: %8.15f\n", dij1i);
	rij = r2*dij1i;
	//printf("rij: %8.15f\n", rij);
	//printf("vdwRad[%d]: %8.15f vdwRad[%d]: %8.15f\n", i, vdwRad[i], j, vdwRad[j]);
	tmpaij = vdwRad[i] - rij*0.5 - (vdwRad[i]*vdwRad[i] - vdwRad[j]*vdwRad[j])*0.5*dij1i;
	//printf("tmpaij: %8.15f\n", tmpaij);
	aij = 2.0*M_PI*vdwRad[i]*tmpaij;
	//printf("aij: %8.8f\n", aij);

	daijddij = M_PI*vdwRad[i]*(dij1i*dij1i*(vdwRad[i]*vdwRad[i] - vdwRad[j]*vdwRad[j]) - 1);
	ddx = -daijddij * xij * dij1i;
	ddy = -daijddij * yij * dij1i;
	ddz = -daijddij * zij * dij1i;

	sumaij = sumaij + aij;
	//printf("sumaij: %8.15f\n", sumaij);

	count2 = icount;
	sumajk2 = 0.0;
	sumkdx = 0.0;
	sumkdy = 0.0;
	sumkdz = 0.0;
	ppa = -surften*(p3[i] + p4[i]*aij);
	//printf("ppa: %8.8f\n", ppa);

	int count2_fin;
	for(int kk = count2; kk < max_count; kk++){
	  if(ineighbor[kk] == -1){
	    count2_fin = kk;
	    break;
	  }
	}

	//printf("count2: %d count2_fin: %d\n", count2, count2_fin);

	num_k_vals = 0;
	for(int kk = count2; kk < count2_fin; kk++){
	  int k = ineighbor[kk];
	  if(k != j){
	    k_vals[num_k_vals] = k;
	    num_k_vals++;
	    //printf("kk: %d k: %d num_k_vals: %d\n", kk, k, num_k_vals);
	  }
	}
	//printf("num_k_vals: %d\n", num_k_vals);

	count2 = icount;

	for(int kk = 0; kk < num_k_vals; kk++){
	  int k = k_vals[kk];
	  //printf("k: %d\n", k);
	  xjk = coords[j].x - coords[k].x;
	  yjk = coords[j].y - coords[k].y;
	  zjk = coords[j].z - coords[k].z;

	  rjk2 = xjk*xjk + yjk*yjk + zjk*zjk;
	  djk1i = 1.0/sqrt(rjk2);
	  rjk = rjk2*djk1i;
	  if((vdwRad[j] + vdwRad[k]) > rjk){
	    vdw2dif = vdwRad[j]*vdwRad[j] - vdwRad[k]*vdwRad[k];
	    tmpajk = 2.0*vdwRad[j] - rjk - vdw2dif*djk1i;
	    ajk = M_PI*vdwRad[j]*tmpajk;
	    //printf("j: %d k: %d ajk: %8.10f\n", j+1, k+1, ajk);
	    sumajk += ajk;
	    //printf("sumajk: %8.10f\n", sumajk);
	    sumajk2 += ajk;

	    dajkddjk = M_PI*vdwRad[j]*djk1i*(djk1i*djk1i*vdw2dif - 1.0);
	    
	    dajx = dajkddjk*xjk;
	    dajy = dajkddjk*yjk;
	    dajz = dajkddjk*zjk;

	    forces[k].x -= dajx*ppa;
	    forces[k].y -= dajy*ppa;
	    forces[k].z -= dajz*ppa;

	    sumkdx += dajx;
	    sumkdy += dajy;
	    sumkdz += dajz;

	    //printf("sumdx: %8.10f\n", sumdx);

	  }
	}

	sumaijajk += aij*sumajk2;
	  
	sumidx -= ddx;
	sumidy -= ddy;
	sumidz -= ddz;

	//printf("ddz: %8.10f sumajk2: %8.10f\n", ddz, sumajk2);
	sumdx -= ddx*sumajk2;
	sumdy -= ddy*sumajk2;
	sumdz -= ddz*sumajk2;

	//printf("sumdz: %8.15f\n", sumdz);

	//printf("ddz: %8.10f sumkdz: %8.10f\n", ddz, sumkdz);
	lastxj = ddx*sumajk2 + aij*sumkdx;
	lastyj = ddy*sumajk2 + aij*sumkdy;
	lastzj = ddz*sumajk2 + aij*sumkdz;

	//printf("lastzj: %8.15f\n", lastzj);

	//printf("ddx: %8.15f sumkdx: %8.15f lastxj: %8.15f\n", ddx, sumkdx, lastxj);
	daidxj = surften*(p2[i]*ddx + p3[i]*sumkdx + p4[i]*lastxj);
	daidyj = surften*(p2[i]*ddy + p3[i]*sumkdy + p4[i]*lastyj);
	daidzj = surften*(p2[i]*ddz + p3[i]*sumkdz + p4[i]*lastzj);

	//printf("daidxj: %8.15f\n", daidxj);
	forces[j].x -= daidxj;
	forces[j].y -= daidyj;
	forces[j].z -= daidzj;
      
	
      }

      //printf("sumaij: %8.10f sumajk: %8.10f sumaijajk: %8.10f\n", sumaij, sumajk, sumaijajk);
      ai = p1[i]*si + p2[i]*sumaij + p3[i]*sumajk + p4[i]*sumaijajk;
      //printf("ai: %8.15f\n", ai);

      daidxi = surften*(p2[i]*sumidx + p4[i]*sumdx);
      daidyi = surften*(p2[i]*sumidy + p4[i]*sumdy);
      daidzi = surften*(p2[i]*sumidz + p4[i]*sumdz);

      forces[i].x -= daidxi;
      forces[i].y -= daidyi;
      forces[i].z -= daidzi;

      totsasa += ai;
      
      //printf("sumidz: %8.10f sumdz: %8.10f daidzi: %8.10f\n", sumidz, sumdz, daidzi);
    }
  }

  esurf = surften*totsasa;
  
  //printf("esurf: %8.15f\n", esurf);

  //printf("egb: %8.10f\n", epol);


  return epol + esurf;
}