rendered paste bodydouble 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;
}