61 #if !defined(VINLINE_VACC) 66 return Vmem_bytes(thee->
mem);
93 VASSERT(thee != VNULL);
98 "Vacc_ivdwAcc: got radius (%g) bigger than max radius (%g)\n",
113 for (iatom=0; iatom<cell->
natoms; iatom++) {
114 atom = cell->
atoms[iatom];
117 if (atom->
id == atomID)
continue;
120 dist2 = VSQR(center[0]-apos[0]) + VSQR(center[1]-apos[1])
121 + VSQR(center[2]-apos[2]);
122 if (dist2 < VSQR(atom->
radius+radius)){
141 thee = (
Vacc*)Vmem_malloc(VNULL, 1,
sizeof(
Vacc) );
142 VASSERT( thee != VNULL);
143 VASSERT(
Vacc_ctor2(thee, alist, clist, surf_density));
161 if (alist == VNULL) {
162 Vnm_print(2,
"Vacc_storeParms: Got NULL Valist!\n");
164 }
else thee->
alist = alist;
165 if (clist == VNULL) {
166 Vnm_print(2,
"Vacc_storeParms: Got NULL Vclist!\n");
168 }
else thee->
clist = clist;
176 if (rad > maxrad) maxrad = rad;
180 maxarea = 4.0*VPI*maxrad*maxrad;
181 nsphere = (int)ceil(maxarea*surf_density);
183 Vnm_print(0,
"Vacc_storeParms: Surf. density = %g\n", surf_density);
184 Vnm_print(0,
"Vacc_storeParms: Max area = %g\n", maxarea);
186 Vnm_print(0,
"Vacc_storeParms: Using %d-point reference sphere\n",
200 thee->
atomFlags = (
int*)Vmem_malloc(thee->
mem, natoms,
sizeof(
int));
203 "Vacc_allocate: Failed to allocate %d (int)s for atomFlags!\n",
207 for (i=0; i<natoms; i++) (thee->
atomFlags)[i] = 0;
221 Vnm_print(2,
"Vacc_ctor2: parameter check failed!\n");
226 thee->
mem = Vmem_ctor(
"APBS::VACC");
227 if (thee->
mem == VNULL) {
228 Vnm_print(2,
"Vacc_ctor2: memory object setup failed!\n");
237 Vnm_print(2,
"Vacc_ctor2: memory allocation failed!\n");
247 if ((*thee) != VNULL) {
249 Vmem_free(VNULL, 1,
sizeof(
Vacc), (
void **)thee);
261 Vmem_free(thee->
mem, natoms,
sizeof(
int), (
void **)&(thee->
atomFlags));
267 if (thee->
surf != VNULL) {
270 (
void **)&(thee->
surf));
274 Vmem_dtor(&(thee->
mem));
292 if (cell == VNULL)
return 1.0;
295 for (iatom=0; iatom<cell->
natoms; iatom++) {
296 atom = cell->
atoms[iatom];
298 dist2 = VSQR(center[0]-apos[0]) + VSQR(center[1]-apos[1])
299 + VSQR(center[2]-apos[2]);
335 VASSERT(thee != NULL);
339 w3i = 1.0/(win*win*win);
342 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
350 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
351 + VSQR(apos[2]-center[2]));
354 if (dist < (arad - win))
return;
357 else if (dist > (arad + win))
return;
361 else if ((VABS(dist - (arad - win)) < VSMALL) ||
362 (VABS(dist - (arad + win)) < VSMALL))
return;
365 sm = dist - arad + win;
367 mychi = 0.75*sm2*w2i -0.25*sm*sm2*w3i;
368 mygrad = 1.5*sm*w2i - 0.75*sm2*w3i;
371 VASSERT(mychi > 0.0);
373 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
396 VASSERT(thee != NULL);
400 w3i = 1.0/(win*win*win);
403 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
411 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
412 + VSQR(apos[2]-center[2]));
415 if (dist < (arad - win))
return;
418 else if (dist > (arad + win))
return;
422 else if ((VABS(dist - (arad - win)) < VSMALL) ||
423 (VABS(dist - (arad + win)) < VSMALL))
return;
426 sm = dist - arad + win;
428 mychi = 0.75*sm2*w2i -0.25*sm*sm2*w3i;
429 mygrad = 1.5*sm*w2i - 0.75*sm2*w3i;
432 VASSERT(mychi > 0.0);
434 grad[i] = -(mygrad)*((center[i] - apos[i])/dist);
456 VASSERT(thee != NULL);
460 w3i = 1.0/(win*win*win);
467 sctot = VMAX2(0, (arad - win));
468 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
469 + VSQR(apos[2]-center[2]));
472 if ((dist < sctot) || (VABS(dist - sctot) < VSMALL)){
475 }
else if ((dist > stot) || (VABS(dist - stot) < VSMALL)) {
479 sm = dist - arad + win;
481 value = 0.75*sm2*w2i - 0.25*sm*sm2*w3i;
506 VASSERT(thee != NULL);
509 for (iatom=0; iatom<cell->
natoms; iatom++) {
511 atom = cell->
atoms[iatom];
520 if (value < VSMALL)
return value;
536 VASSERT(thee != NULL);
539 Vnm_print(2,
"Vacc_splineAcc: Vclist has max_radius=%g;\n",
541 Vnm_print(2,
"Vacc_splineAcc: Insufficient for win=%g, infrad=%g\n",
548 if (cell == VNULL)
return 1.0;
552 for (iatom=0; iatom<cell->
natoms; iatom++) {
553 atom = cell->
atoms[iatom];
558 return splineAcc(thee, center, win, infrad, cell);
562 double win,
double infrad,
double *grad) {
564 int iatom, i, atomID;
570 VASSERT(thee != NULL);
573 Vnm_print(2,
"Vacc_splineAccGrad: Vclist max_radius=%g;\n",
575 Vnm_print(2,
"Vacc_splineAccGrad: Insufficient for win=%g, infrad=%g\n",
581 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
585 if (cell == VNULL)
return;
588 for (iatom=0; iatom<cell->
natoms; iatom++) {
589 atom = cell->
atoms[iatom];
595 acc =
splineAcc(thee, center, win, infrad, cell);
599 for (iatom=0; iatom<cell->
natoms; iatom++) {
600 atom = cell->
atoms[iatom];
603 for (i=0; i<
VAPBS_DIM; i++) grad[i] += tgrad[i];
605 for (i=0; i<
VAPBS_DIM; i++) grad[i] *= -acc;
643 int ipt, iatom, atomID;
646 rad2 = radius*radius;
654 Vnm_print(2,
"Vacc_fastMolAcc: unexpected VNULL VclistCell!\n");
659 for (iatom=0; iatom<cell->
natoms; iatom++) {
660 atom = cell->
atoms[iatom];
662 surf = thee->
surf[atomID];
664 for (ipt=0; ipt<surf->
npts; ipt++) {
666 dist2 = VSQR(center[0]-(surf->
xpts[ipt]))
667 + VSQR(center[1]-(surf->
ypts[ipt]))
668 + VSQR(center[2]-(surf->
zpts[ipt]));
669 if (dist2 < rad2)
return 1.0;
678 #if defined(HAVE_MC_H) 679 VPUBLIC
void Vacc_writeGMV(
Vacc *thee,
double radius,
int meth, Gem *gm,
680 char *iodev,
char *iofmt,
char *iohost,
char *iofile) {
682 double *accVals[MAXV], coord[3];
686 for (ivert=0; ivert<MAXV; ivert++) accVals[ivert] = VNULL;
687 accVals[0] = (
void *)Vmem_malloc(thee->
mem, Gem_numVV(gm),
sizeof(double));
688 accVals[1] = (
void *)Vmem_malloc(thee->
mem, Gem_numVV(gm),
sizeof(double));
689 for (ivert=0; ivert<Gem_numVV(gm); ivert++) {
690 for (icoord=0;icoord<3;icoord++)
691 coord[icoord] = VV_coord(Gem_VV(gm, ivert), icoord);
693 accVals[0][ivert] =
Vacc_molAcc(thee, coord, radius);
694 accVals[1][ivert] =
Vacc_molAcc(thee, coord, radius);
695 }
else if (meth == 1) {
698 }
else if (meth == 2) {
703 sock = Vio_ctor(iodev, iofmt, iohost, iofile,
"w");
704 Gem_writeGMV(gm, sock, 1, accVals);
706 Vmem_free(thee->
mem, Gem_numVV(gm),
sizeof(
double),
707 (
void **)&(accVals[0]));
708 Vmem_free(thee->
mem, Gem_numVV(gm),
sizeof(
double),
709 (
void **)&(accVals[1]));
732 if (thee->
surf == VNULL) {
735 #if defined(DEBUG_MAC_OSX_OCL) || defined(DEBUG_MAC_OSX_STANDARD) 736 #include "mach_chud.h" 738 #pragma omp parallel for private(i,atom) 740 for (i=0; i<natom; i++) {
751 for (i=0; i<natom; i++) {
753 asurf = thee->
surf[i];
756 Vnm_print(2,
"Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
760 asurf = thee->
surf[i];
762 area += (asurf->
area);
765 #if defined(DEBUG_MAC_OSX_OCL) || defined(DEBUG_MAC_OSX_STANDARD) 766 mets_(&mbeg,
"Vacc_SASA - Parallel");
769 Vnm_print(0,
"Vacc_SASA: Time elapsed: %f\n", ((
double)clock() - ts) / CLOCKS_PER_SEC);
788 asurf = thee->
surf[id];
792 Vnm_print(2,
"Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
796 asurf = thee->
surf[id];
820 thee->
npts = nsphere;
824 if (thee->
npts > 0) {
825 thee->
xpts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
double));
826 thee->
ypts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
double));
827 thee->
zpts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
double));
828 thee->
bpts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
char));
843 if ((*thee) != VNULL) {
855 if (thee->
npts > 0) {
856 Vmem_free(thee->
mem, thee->
npts,
sizeof(
double),
857 (
void **)&(thee->
xpts));
858 Vmem_free(thee->
mem, thee->
npts,
sizeof(
double),
859 (
void **)&(thee->
ypts));
860 Vmem_free(thee->
mem, thee->
npts,
sizeof(
double),
861 (
void **)&(thee->
zpts));
862 Vmem_free(thee->
mem, thee->
npts,
sizeof(
char),
863 (
void **)&(thee->
bpts));
874 double arad, rad, pos[3], *apos;
889 for (i=0; i<ref->
npts; i++) {
891 pos[0] = rad*(ref->
xpts[i]) + apos[0];
892 pos[1] = rad*(ref->
ypts[i]) + apos[1];
893 pos[2] = rad*(ref->
zpts[i]) + apos[2];
896 ref->
bpts[i] = (ref->
bpts[i] << 1) + 1;
907 for (i=0; i<ref->
npts; i++) {
908 char flag = ref->
bpts[i] & 1;
912 surf->
xpts[j] = rad*(ref->
xpts[i]) + apos[0];
913 surf->
ypts[j] = rad*(ref->
ypts[i]) + apos[1];
914 surf->
zpts[j] = rad*(ref->
zpts[i]) + apos[2];
920 surf->
area = 4.0*VPI*rad*rad*((double)(surf->
npts))/((
double)(ref->
npts));
929 int nactual, i, itheta, ntheta, iphi, nphimax, nphi;
931 double sintheta, costheta, theta, dtheta;
932 double sinphi, cosphi, phi, dphi;
935 frac = ((double)(npts))/4.0;
936 ntheta = VRINT(VSQRT(
Vunit_pi*frac));
937 dtheta =
Vunit_pi/((double)(ntheta));
942 for (itheta=0; itheta<ntheta; itheta++) {
943 theta = dtheta*((double)(itheta));
944 sintheta = VSIN(theta);
945 costheta = VCOS(theta);
946 nphi = VRINT(sintheta*nphimax);
954 for (i=0; i<nactual; i++) surf->
bpts[i] = 1;
958 for (itheta=0; itheta<ntheta; itheta++) {
959 theta = dtheta*((double)(itheta));
960 sintheta = VSIN(theta);
961 costheta = VCOS(theta);
962 nphi = VRINT(sintheta*nphimax);
965 for (iphi=0; iphi<nphi; iphi++) {
966 phi = dphi*((double)(iphi));
969 surf->
xpts[nactual] = cosphi * sintheta;
970 surf->
ypts[nactual] = sinphi * sintheta;
971 surf->
zpts[nactual] = costheta;
977 surf->
npts = nactual;
991 asurf = thee->
surf[id];
995 Vnm_print(2,
"Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
999 asurf = thee->
surf[id];
1007 double win,
double infrad,
Vatom *atom,
double *grad) {
1010 double dist, *apos, arad, sm, sm2, sm3, sm4, sm5, sm6, sm7;
1011 double e, e2, e3, e4, e5, e6, e7;
1012 double b, b2, b3, b4, b5, b6, b7;
1013 double c0, c1, c2, c3, c4, c5, c6, c7;
1014 double denom, mygrad;
1017 VASSERT(thee != NULL);
1020 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
1029 arad = arad + infrad;
1046 denom = e7 - 7.0*b*e6 + 21.0*b2*e5 - 35.0*e4*b3
1047 + 35.0*e3*b4 - 21.0*b5*e2 + 7.0*e*b6 - b7;
1048 c0 = b4*(35.0*e3 - 21.0*b*e2 + 7*e*b2 - b3)/denom;
1049 c1 = -140.0*b3*e3/denom;
1050 c2 = 210.0*e2*b2*(e + b)/denom;
1051 c3 = -140.0*e*b*(e2 + 3.0*b*e + b2)/denom;
1052 c4 = 35.0*(e3 + 9.0*b*e2 + + 9.0*e*b2 + b3)/denom;
1053 c5 = -84.0*(e2 + 3.0*b*e + b2)/denom;
1054 c6 = 70.0*(e + b)/denom;
1057 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1058 + VSQR(apos[2]-center[2]));
1062 if (dist < (arad - win))
return;
1065 else if (dist > (arad + win))
return;
1069 else if ((VABS(dist - (arad - win)) < VSMALL) ||
1070 (VABS(dist - (arad + win)) < VSMALL))
return;
1080 mychi = c0 + c1*sm + c2*sm2 + c3*sm3
1081 + c4*sm4 + c5*sm5 + c6*sm6 + c7*sm7;
1082 mygrad = c1 + 2.0*c2*sm + 3.0*c3*sm2 + 4.0*c4*sm3
1083 + 5.0*c5*sm4 + 6.0*c6*sm5 + 7.0*c7*sm6;
1087 }
else if (mychi > 1.0) {
1093 VASSERT(mychi > 0.0);
1095 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1100 double win,
double infrad,
Vatom *atom,
double *grad) {
1103 double dist, *apos, arad, sm, sm2, sm3, sm4, sm5;
1104 double e, e2, e3, e4, e5;
1105 double b, b2, b3, b4, b5;
1106 double c0, c1, c2, c3, c4, c5;
1107 double denom, mygrad;
1110 VASSERT(thee != NULL);
1113 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
1122 arad = arad + infrad;
1135 denom = pow((e - b), 5.0);
1136 c0 = -10.0*e2*b3 + 5.0*e*b4 - b5;
1138 c2 = -30.0*(e2*b + e*b2);
1139 c3 = 10.0*(e2 + 4.0*e*b + b2);
1149 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1150 + VSQR(apos[2]-center[2]));
1154 if (dist < (arad - win))
return;
1157 else if (dist > (arad + win))
return;
1161 else if ((VABS(dist - (arad - win)) < VSMALL) ||
1162 (VABS(dist - (arad + win)) < VSMALL))
return;
1170 mychi = c0 + c1*sm + c2*sm2 + c3*sm3
1172 mygrad = c1 + 2.0*c2*sm + 3.0*c3*sm2 + 4.0*c4*sm3
1177 }
else if (mychi > 1.0) {
1183 VASSERT(mychi > 0.0);
1185 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1226 if(tRad == 0.0)
return;
1228 area = 4.0*VPI*(tRad+srad)*(tRad+srad)/((double)(ref->npts));
1229 for (ipt=0; ipt<ref->npts; ipt++) {
1230 vec[0] = (tRad+srad)*ref->xpts[ipt] + tPos[0];
1231 vec[1] = (tRad+srad)*ref->ypts[ipt] + tPos[1];
1232 vec[2] = (tRad+srad)*ref->zpts[ipt] + tPos[2];
1234 dx = dx+vec[0]-tPos[0];
1235 dy = dy+vec[1]-tPos[1];
1236 dz = dz+vec[2]-tPos[2];
1240 if ((tRad+srad) != 0){
1241 dSA[0] = dx*area/(tRad+srad);
1242 dSA[1] = dy*area/(tRad+srad);
1243 dSA[2] = dz*area/(tRad+srad);
1252 VPRIVATE
double Vacc_SASAPos(
Vacc *thee,
double radius) {
1263 for (i=0; i<natom; i++) {
1265 asurf = thee->
surf[i];
1269 asurf = thee->
surf[i];
1270 area += (asurf->
area);
1277 VPRIVATE
double Vacc_atomSASAPos(
Vacc *thee,
1285 static int warned = 0;
1287 if ((thee->
surf == VNULL) || (mode == 1)){
1289 Vnm_print(2,
"WARNING: Recalculating entire surface!!!!\n");
1292 Vacc_SASAPos(thee, radius);
1296 asurf = thee->
surf[id];
1300 asurf = thee->
surf[id];
1345 tPos[0] = temp_Pos[0];
1346 tPos[1] = temp_Pos[1];
1347 tPos[2] = temp_Pos[2];
1350 temp_Pos[0] -= dpos;
1351 axb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1352 temp_Pos[0] = tPos[0];
1354 temp_Pos[0] += dpos;
1355 axt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1356 temp_Pos[0] = tPos[0];
1359 temp_Pos[1] -= dpos;
1360 ayb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1361 temp_Pos[1] = tPos[1];
1363 temp_Pos[1] += dpos;
1364 ayt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1365 temp_Pos[1] = tPos[1];
1368 temp_Pos[2] -= dpos;
1369 azb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1370 temp_Pos[2] = tPos[2];
1372 temp_Pos[2] += dpos;
1373 azt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1374 temp_Pos[2] = tPos[2];
1377 Vacc_atomSASAPos(thee, srad, atom,0);
1380 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1381 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1382 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1392 double *temp_Pos, tRad;
1394 double axb1,axt1,ayb1,ayt1,azb1,azt1;
1407 tPos[0] = temp_Pos[0];
1408 tPos[1] = temp_Pos[1];
1409 tPos[2] = temp_Pos[2];
1412 temp_Pos[0] -= dpos;
1413 axb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1414 temp_Pos[0] = tPos[0];
1416 temp_Pos[0] += dpos;
1417 axt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1418 temp_Pos[0] = tPos[0];
1421 temp_Pos[1] -= dpos;
1422 ayb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1423 temp_Pos[1] = tPos[1];
1425 temp_Pos[1] += dpos;
1426 ayt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1427 temp_Pos[1] = tPos[1];
1430 temp_Pos[2] -= dpos;
1431 azb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1432 temp_Pos[2] = tPos[2];
1434 temp_Pos[2] += dpos;
1435 azt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1436 temp_Pos[2] = tPos[2];
1439 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1440 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1441 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1451 double *temp_Pos, tRad;
1453 double axb1,axt1,ayb1,ayt1,azb1,azt1;
1466 tPos[0] = temp_Pos[0];
1467 tPos[1] = temp_Pos[1];
1468 tPos[2] = temp_Pos[2];
1471 temp_Pos[0] -= dpos;
1473 temp_Pos[0] = tPos[0];
1475 temp_Pos[0] += dpos;
1477 temp_Pos[0] = tPos[0];
1480 temp_Pos[1] -= dpos;
1482 temp_Pos[1] = tPos[1];
1484 temp_Pos[1] += dpos;
1486 temp_Pos[1] = tPos[1];
1489 temp_Pos[2] -= dpos;
1491 temp_Pos[2] = tPos[2];
1493 temp_Pos[2] += dpos;
1495 temp_Pos[2] = tPos[2];
1498 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1499 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1500 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1508 double spacs[3], vec[3];
1509 double w, wx, wy, wz, len, fn, x, y, z, vol;
1510 double vol_density,sav;
1511 double *lower_corner, *upper_corner;
1520 for (i=0; i<3; i++) {
1521 len = upper_corner[i] - lower_corner[i];
1523 fn = len*vol_density + 1;
1524 npts[i] = (int)ceil(fn);
1525 spacs[i] = len/((double)(npts[i])-1.0);
1526 if (apolparm != VNULL) {
1528 if (apolparm->
grid[i] > spacs[i]) {
1529 Vnm_print(2,
"Vacc_totalSAV: Warning, your GRID value (%g) is larger than the recommended value (%g)!\n",
1530 apolparm->
grid[i], spacs[i]);
1532 spacs[i] = apolparm->
grid[i];
1538 for (x=lower_corner[0]; x<=upper_corner[0]; x=x+spacs[0]) {
1539 if ( VABS(x - lower_corner[0]) < VSMALL) {
1541 }
else if ( VABS(x - upper_corner[0]) < VSMALL) {
1547 for (y=lower_corner[1]; y<=upper_corner[1]; y=y+spacs[1]) {
1548 if ( VABS(y - lower_corner[1]) < VSMALL) {
1550 }
else if ( VABS(y - upper_corner[1]) < VSMALL) {
1556 for (z=lower_corner[2]; z<=upper_corner[2]; z=z+spacs[2]) {
1557 if ( VABS(z - lower_corner[2]) < VSMALL) {
1559 }
else if ( VABS(z - upper_corner[2]) < VSMALL) {
1574 w = spacs[0]*spacs[1]*spacs[2];
1581 Vclist *clist,
int iatom,
double *value) {
1587 int xmin, ymin, zmin;
1588 int xmax, ymax, zmax;
1590 double sigma6, sigma12;
1592 double spacs[3], vec[3];
1593 double w, wx, wy, wz, len, fn, x, y, z, vol;
1595 double vol_density, energy, rho, srad;
1596 double psig, epsilon, watepsilon, sigma, watsigma, eni, chi;
1599 double *lower_corner, *upper_corner;
1601 Vatom *atom = VNULL;
1602 VASSERT(apolparm != VNULL);
1621 srad = apolparm->
srad;
1622 rho = apolparm->
bconc;
1627 sigma = psig + watsigma;
1628 epsilon = VSQRT((epsilon * watepsilon));
1631 sigma6 = VPOW(sigma,6);
1632 sigma12 = VPOW(sigma,12);
1635 xmin = pos[0] - pad;
1636 xmax = pos[0] + pad;
1637 ymin = pos[1] - pad;
1638 ymax = pos[1] + pad;
1639 zmin = pos[2] - pad;
1640 zmax = pos[2] + pad;
1642 for (i=0; i<3; i++) {
1643 len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1645 fn = len*vol_density + 1;
1646 npts[i] = (int)ceil(fn);
1649 if (apolparm->
grid[i] > spacs[i]) {
1650 Vnm_print(2,
"Vacc_totalSAV: Warning, your GRID value (%g) is larger than the recommended value (%g)!\n",
1651 apolparm->
grid[i], spacs[i]);
1653 spacs[i] = apolparm->
grid[i];
1657 for (x=xmin; x<=xmax; x=x+spacs[0]) {
1658 if ( VABS(x - xmin) < VSMALL) {
1660 }
else if ( VABS(x - xmax) < VSMALL) {
1666 for (y=ymin; y<=ymax; y=y+spacs[1]) {
1667 if ( VABS(y - ymin) < VSMALL) {
1669 }
else if ( VABS(y - ymax) < VSMALL) {
1675 for (z=zmin; z<=zmax; z=z+spacs[2]) {
1676 if ( VABS(z - zmin) < VSMALL) {
1678 }
else if ( VABS(z - zmax) < VSMALL) {
1689 if (VABS(chi) > VSMALL) {
1691 x2 = VSQR(vec[0]-pos[0]);
1692 y2 = VSQR(vec[1]-pos[1]);
1693 z2 = VSQR(vec[2]-pos[2]);
1694 r = VSQRT(x2+y2+z2);
1696 if (r <= 14 && r >= sigma) {
1697 eni = chi*rho*epsilon*(-2.0*sigma6/VPOW(r,6)+sigma12/VPOW(r,12));
1699 eni = -1.0*epsilon*chi*rho;
1713 w = spacs[0]*spacs[1]*spacs[2];
1727 double energy = 0.0;
1728 double tenergy = 0.0;
1729 double rho = apolparm->
bconc;
1733 if(apolparm->
setwat == 0){
1734 Vnm_print(2,
"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1738 if (VABS(rho) < VSMALL) {
1745 if(rc == 0)
return 0;
1810 VASSERT(apolparm != VNULL);
1814 if(apolparm->
setwat == 0){
1815 Vnm_print(2,
"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1827 srad = apolparm->
srad;
1828 rho = apolparm->
bconc;
1834 sigma = psig + watsigma;
1835 epsilon = VSQRT((epsilon * watepsilon));
1838 sigma6 = VPOW(sigma,6);
1839 sigma12 = VPOW(sigma,12);
1842 for (i=0; i<3; i++) {
1843 len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1845 fn = len*vol_density + 1;
1846 npts[i] = (int)ceil(fn);
1850 if (apolparm->
grid[i] > spacs[i]) {
1851 Vnm_print(2,
"Vacc_totalSAV: Warning, your GRID value (%g) is larger than the recommended value (%g)!\n",
1852 apolparm->
grid[i], spacs[i]);
1854 spacs[i] = apolparm->
grid[i];
1858 xmin = pos[0] - pad;
1859 xmax = pos[0] + pad;
1860 ymin = pos[1] - pad;
1861 ymax = pos[1] + pad;
1862 zmin = pos[2] - pad;
1863 zmax = pos[2] + pad;
1865 for (x=xmin; x<=xmax; x=x+spacs[0]) {
1866 if ( VABS(x - xmin) < VSMALL) {
1868 }
else if ( VABS(x - xmax) < VSMALL) {
1874 for (y=ymin; y<=ymax; y=y+spacs[1]) {
1875 if ( VABS(y - ymin) < VSMALL) {
1877 }
else if ( VABS(y - ymax) < VSMALL) {
1883 for (z=zmin; z<=zmax; z=z+spacs[2]) {
1884 if ( VABS(z - zmin) < VSMALL) {
1886 }
else if ( VABS(z - zmax) < VSMALL) {
1899 x2 = VSQR(vec[0]-pos[0]);
1900 y2 = VSQR(vec[1]-pos[1]);
1901 z2 = VSQR(vec[2]-pos[2]);
1902 r = VSQRT(x2+y2+z2);
1904 if (r <= 14 && r >= sigma){
1906 fo = 12.0*chi*rho*epsilon*(sigma6/VPOW(r,7)-sigma12/VPOW(r,13));
1908 fpt[0] = -1.0*(pos[0]-vec[0])*fo/r;
1909 fpt[1] = -1.0*(pos[1]-vec[1])*fo/r;
1910 fpt[2] = -1.0*(pos[2]-vec[2])*fo/r;
1913 for (si=0; si < 3; si++) fpt[si] = 0.0;
1916 for (si=0; si < 3; si++) fpt[si] = 0.0;
1920 force[i] += (w*fpt[i]);
1927 w = spacs[0]*spacs[1]*spacs[2];
1928 for(i=0;i<3;i++) force[i] *= w;
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist, Vclist *clist, int iatom, double *value)
Calculate the WCA energy for an atom.
VPUBLIC void Vacc_totalAtomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Testing purposes only.
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
double upper_corner[VAPBS_DIM]
Oracle for solvent- and ion-accessibility around a biomolecule.
VPRIVATE int Vacc_storeParms(Vacc *thee, Valist *alist, Vclist *clist, double surf_density)
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
VPRIVATE int ivdwAccExclus(Vacc *thee, double center[3], double radius, int atomID)
Determines if a point is within the union of the spheres centered at the atomic centers with radii eq...
VPUBLIC double Vclist_maxRadius(Vclist *thee)
Get the max probe radius value (in A) the cell list was constructed with.
Contains declarations for class Vacc.
VPUBLIC unsigned long int Vacc_memChk(Vacc *thee)
Get number of bytes in this object and its members.
VPUBLIC void Vacc_splineAccGradAtomNorm(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, Vatom *atom, double *grad)
Report gradient of spline-based accessibility with respect to a particular atom normalized by the acc...
VEXTERNC double Vacc_ivdwAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report inflated van der Waals accessibility.
VPUBLIC void Vacc_splineAccGradAtomUnnorm(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, Vatom *atom, double *grad)
Report gradient of spline-based accessibility with respect to a particular atom (see Vpmg_splineAccAt...
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
VPUBLIC void Vacc_splineAccGradAtomNorm4(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, Vatom *atom, double *grad)
Report gradient of spline-based accessibility with respect to a particular atom normalized by a 4th o...
VPUBLIC void VaccSurf_dtor2(VaccSurf *thee)
Destroy the surface object.
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
VPUBLIC VaccSurf * VaccSurf_refSphere(Vmem *mem, int npts)
Set up an array of points for a reference sphere of unit radius.
VPUBLIC VclistCell * Vclist_getCell(Vclist *thee, double pos[VAPBS_DIM])
Return cell corresponding to specified position or return VNULL.
VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
VPUBLIC int Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
VPUBLIC void Vacc_dtor2(Vacc *thee)
FORTRAN stub to destroy object.
VPUBLIC void Vacc_splineAccGradAtomNorm3(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, Vatom *atom, double *grad)
Report gradient of spline-based accessibility with respect to a particular atom normalized by a 3rd o...
VPUBLIC VaccSurf * VaccSurf_ctor(Vmem *mem, double probe_radius, int nsphere)
Allocate and construct the surface object; do not assign surface points to positions.
VEXTERNC double Vacc_vdwAcc(Vacc *thee, double center[VAPBS_DIM])
Report van der Waals accessibility.
VPUBLIC void Vacc_splineAccGrad(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, double *grad)
Report gradient of spline-based accessibility.
VPUBLIC double Vacc_splineAccAtom(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, Vatom *atom)
Report spline-based accessibility for a given atom.
VPRIVATE double splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, VclistCell *cell)
Fast spline-based surface computation subroutine.
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
VPUBLIC double Vatom_getAtomID(Vatom *thee)
Get atom ID.
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
VPUBLIC void Vacc_totalAtomdSAV(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA, Vclist *clist)
Total solvent accessible volume.
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
VPUBLIC double Vacc_fastMolAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility quickly.
VPUBLIC double Vacc_SASA(Vacc *thee, double radius)
Build the solvent accessible surface (SAS) and calculate the solvent accessible surface area.
VPUBLIC int VaccSurf_ctor2(VaccSurf *thee, Vmem *mem, double probe_radius, int nsphere)
Construct the surface object using previously allocated memory; do not assign surface points to posit...
#define VAPBS_DIM
Our dimension.
VPUBLIC VaccSurf * Vacc_atomSASPoints(Vacc *thee, double radius, Vatom *atom)
Get the set of points for this atom's solvent-accessible surface.
Surface object list of per-atom surface points.
VPUBLIC void VaccSurf_dtor(VaccSurf **thee)
Destroy the surface object and free its memory.
VPRIVATE int Vacc_allocate(Vacc *thee)
VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility.
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
VPUBLIC int Vacc_ctor2(Vacc *thee, Valist *alist, Vclist *clist, double surf_density)
FORTRAN stub to construct the accessibility object.
Contains public data members for Vatom class/module.
Container class for list of atom objects.
double lower_corner[VAPBS_DIM]
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad)
Report spline-based accessibility.
Parameter structure for APOL-specific variables from input files.
VPUBLIC VaccSurf * Vacc_atomSurf(Vacc *thee, Vatom *atom, VaccSurf *ref, double prad)
Set up an array of points corresponding to the SAS due to a particular atom.