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];
808 Vnm_print(2,
"VaccSurf_ctor: Error! The requested number of grid points (%d) exceeds the maximum (%d)!\n", nsphere,
MAX_SPHERE_PTS);
809 Vnm_print(2,
"VaccSurf_ctor: Please check the variable MAX_SPHERE_PTS to reset.\n");
825 thee->
npts = nsphere;
829 if (thee->
npts > 0) {
830 thee->
xpts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
double));
831 thee->
ypts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
double));
832 thee->
zpts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
double));
833 thee->
bpts = Vmem_malloc(thee->
mem, thee->
npts,
sizeof(
char));
848 if ((*thee) != VNULL) {
860 if (thee->
npts > 0) {
861 Vmem_free(thee->
mem, thee->
npts,
sizeof(
double),
862 (
void **)&(thee->
xpts));
863 Vmem_free(thee->
mem, thee->
npts,
sizeof(
double),
864 (
void **)&(thee->
ypts));
865 Vmem_free(thee->
mem, thee->
npts,
sizeof(
double),
866 (
void **)&(thee->
zpts));
867 Vmem_free(thee->
mem, thee->
npts,
sizeof(
char),
868 (
void **)&(thee->
bpts));
903 for (i=0; i<ref->
npts; i++) {
905 pos[0] = rad*(ref->
xpts[i]) + apos[0];
906 pos[1] = rad*(ref->
ypts[i]) + apos[1];
907 pos[2] = rad*(ref->
zpts[i]) + apos[2];
921 for (i=0; i<ref->
npts; i++) {
924 surf->
xpts[j] = rad*(ref->
xpts[i]) + apos[0];
925 surf->
ypts[j] = rad*(ref->
ypts[i]) + apos[1];
926 surf->
zpts[j] = rad*(ref->
zpts[i]) + apos[2];
932 surf->
area = 4.0*VPI*rad*rad*((double)(surf->
npts))/((
double)(ref->
npts));
941 int nactual, i, itheta, ntheta, iphi, nphimax, nphi;
943 double sintheta, costheta, theta, dtheta;
944 double sinphi, cosphi, phi, dphi;
947 frac = ((double)(npts))/4.0;
948 ntheta = VRINT(VSQRT(
Vunit_pi*frac));
949 dtheta =
Vunit_pi/((double)(ntheta));
954 for (itheta=0; itheta<ntheta; itheta++) {
955 theta = dtheta*((double)(itheta));
956 sintheta = VSIN(theta);
957 costheta = VCOS(theta);
958 nphi = VRINT(sintheta*nphimax);
966 for (i=0; i<nactual; i++) surf->
bpts[i] = 1;
970 for (itheta=0; itheta<ntheta; itheta++) {
971 theta = dtheta*((double)(itheta));
972 sintheta = VSIN(theta);
973 costheta = VCOS(theta);
974 nphi = VRINT(sintheta*nphimax);
977 for (iphi=0; iphi<nphi; iphi++) {
978 phi = dphi*((double)(iphi));
981 surf->
xpts[nactual] = cosphi * sintheta;
982 surf->
ypts[nactual] = sinphi * sintheta;
983 surf->
zpts[nactual] = costheta;
989 surf->
npts = nactual;
1003 asurf = thee->
surf[id];
1007 Vnm_print(2,
"Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
1011 asurf = thee->
surf[id];
1019 double win,
double infrad,
Vatom *atom,
double *grad) {
1022 double dist, *apos, arad, sm, sm2, sm3, sm4, sm5, sm6, sm7;
1023 double e, e2, e3, e4, e5, e6, e7;
1024 double b, b2, b3, b4, b5, b6, b7;
1025 double c0, c1, c2, c3, c4, c5, c6, c7;
1026 double denom, mygrad;
1029 VASSERT(thee != NULL);
1032 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
1041 arad = arad + infrad;
1058 denom = e7 - 7.0*b*e6 + 21.0*b2*e5 - 35.0*e4*b3
1059 + 35.0*e3*b4 - 21.0*b5*e2 + 7.0*e*b6 - b7;
1060 c0 = b4*(35.0*e3 - 21.0*b*e2 + 7*e*b2 - b3)/denom;
1061 c1 = -140.0*b3*e3/denom;
1062 c2 = 210.0*e2*b2*(e + b)/denom;
1063 c3 = -140.0*e*b*(e2 + 3.0*b*e + b2)/denom;
1064 c4 = 35.0*(e3 + 9.0*b*e2 + + 9.0*e*b2 + b3)/denom;
1065 c5 = -84.0*(e2 + 3.0*b*e + b2)/denom;
1066 c6 = 70.0*(e + b)/denom;
1069 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1070 + VSQR(apos[2]-center[2]));
1074 if (dist < (arad - win))
return;
1077 else if (dist > (arad + win))
return;
1081 else if ((VABS(dist - (arad - win)) < VSMALL) ||
1082 (VABS(dist - (arad + win)) < VSMALL))
return;
1092 mychi = c0 + c1*sm + c2*sm2 + c3*sm3
1093 + c4*sm4 + c5*sm5 + c6*sm6 + c7*sm7;
1094 mygrad = c1 + 2.0*c2*sm + 3.0*c3*sm2 + 4.0*c4*sm3
1095 + 5.0*c5*sm4 + 6.0*c6*sm5 + 7.0*c7*sm6;
1099 }
else if (mychi > 1.0) {
1105 VASSERT(mychi > 0.0);
1107 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1112 double win,
double infrad,
Vatom *atom,
double *grad) {
1115 double dist, *apos, arad, sm, sm2, sm3, sm4, sm5;
1116 double e, e2, e3, e4, e5;
1117 double b, b2, b3, b4, b5;
1118 double c0, c1, c2, c3, c4, c5;
1119 double denom, mygrad;
1122 VASSERT(thee != NULL);
1125 for (i=0; i<
VAPBS_DIM; i++) grad[i] = 0.0;
1134 arad = arad + infrad;
1147 denom = pow((e - b), 5.0);
1148 c0 = -10.0*e2*b3 + 5.0*e*b4 - b5;
1150 c2 = -30.0*(e2*b + e*b2);
1151 c3 = 10.0*(e2 + 4.0*e*b + b2);
1161 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1162 + VSQR(apos[2]-center[2]));
1166 if (dist < (arad - win))
return;
1169 else if (dist > (arad + win))
return;
1173 else if ((VABS(dist - (arad - win)) < VSMALL) ||
1174 (VABS(dist - (arad + win)) < VSMALL))
return;
1182 mychi = c0 + c1*sm + c2*sm2 + c3*sm3
1184 mygrad = c1 + 2.0*c2*sm + 3.0*c3*sm2 + 4.0*c4*sm3
1189 }
else if (mychi > 1.0) {
1195 VASSERT(mychi > 0.0);
1197 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1238 if(tRad == 0.0)
return;
1240 area = 4.0*VPI*(tRad+
srad)*(tRad+
srad)/((double)(ref->npts));
1241 for (ipt=0; ipt<ref->npts; ipt++) {
1242 vec[0] = (tRad+
srad)*ref->xpts[ipt] + tPos[0];
1243 vec[1] = (tRad+
srad)*ref->ypts[ipt] + tPos[1];
1244 vec[2] = (tRad+
srad)*ref->zpts[ipt] + tPos[2];
1246 dx = dx+vec[0]-tPos[0];
1247 dy = dy+vec[1]-tPos[1];
1248 dz = dz+vec[2]-tPos[2];
1252 if ((tRad+
srad) != 0){
1253 dSA[0] = dx*area/(tRad+
srad);
1254 dSA[1] = dy*area/(tRad+
srad);
1255 dSA[2] = dz*area/(tRad+
srad);
1264 VPRIVATE
double Vacc_SASAPos(
Vacc *thee,
double radius) {
1275 for (i=0; i<natom; i++) {
1277 asurf = thee->
surf[i];
1281 asurf = thee->
surf[i];
1282 area += (asurf->
area);
1289 VPRIVATE
double Vacc_atomSASAPos(
Vacc *thee,
1297 static int warned = 0;
1299 if ((thee->
surf == VNULL) || (mode == 1)){
1301 Vnm_print(2,
"WARNING: Recalculating entire surface!!!!\n");
1304 Vacc_SASAPos(thee, radius);
1308 asurf = thee->
surf[id];
1312 asurf = thee->
surf[id];
1357 tPos[0] = temp_Pos[0];
1358 tPos[1] = temp_Pos[1];
1359 tPos[2] = temp_Pos[2];
1362 temp_Pos[0] -= dpos;
1363 axb1 = Vacc_atomSASAPos(thee,
srad, atom,0);
1364 temp_Pos[0] = tPos[0];
1366 temp_Pos[0] += dpos;
1367 axt1 = Vacc_atomSASAPos(thee,
srad, atom,0);
1368 temp_Pos[0] = tPos[0];
1371 temp_Pos[1] -= dpos;
1372 ayb1 = Vacc_atomSASAPos(thee,
srad, atom,0);
1373 temp_Pos[1] = tPos[1];
1375 temp_Pos[1] += dpos;
1376 ayt1 = Vacc_atomSASAPos(thee,
srad, atom,0);
1377 temp_Pos[1] = tPos[1];
1380 temp_Pos[2] -= dpos;
1381 azb1 = Vacc_atomSASAPos(thee,
srad, atom,0);
1382 temp_Pos[2] = tPos[2];
1384 temp_Pos[2] += dpos;
1385 azt1 = Vacc_atomSASAPos(thee,
srad, atom,0);
1386 temp_Pos[2] = tPos[2];
1389 Vacc_atomSASAPos(thee,
srad, atom,0);
1392 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1393 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1394 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1404 double *temp_Pos, tRad;
1406 double axb1,axt1,ayb1,ayt1,azb1,azt1;
1419 tPos[0] = temp_Pos[0];
1420 tPos[1] = temp_Pos[1];
1421 tPos[2] = temp_Pos[2];
1424 temp_Pos[0] -= dpos;
1425 axb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1426 temp_Pos[0] = tPos[0];
1428 temp_Pos[0] += dpos;
1429 axt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1430 temp_Pos[0] = tPos[0];
1433 temp_Pos[1] -= dpos;
1434 ayb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1435 temp_Pos[1] = tPos[1];
1437 temp_Pos[1] += dpos;
1438 ayt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1439 temp_Pos[1] = tPos[1];
1442 temp_Pos[2] -= dpos;
1443 azb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1444 temp_Pos[2] = tPos[2];
1446 temp_Pos[2] += dpos;
1447 azt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1448 temp_Pos[2] = tPos[2];
1451 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1452 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1453 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1463 double *temp_Pos, tRad;
1465 double axb1,axt1,ayb1,ayt1,azb1,azt1;
1478 tPos[0] = temp_Pos[0];
1479 tPos[1] = temp_Pos[1];
1480 tPos[2] = temp_Pos[2];
1483 temp_Pos[0] -= dpos;
1485 temp_Pos[0] = tPos[0];
1487 temp_Pos[0] += dpos;
1489 temp_Pos[0] = tPos[0];
1492 temp_Pos[1] -= dpos;
1494 temp_Pos[1] = tPos[1];
1496 temp_Pos[1] += dpos;
1498 temp_Pos[1] = tPos[1];
1501 temp_Pos[2] -= dpos;
1503 temp_Pos[2] = tPos[2];
1505 temp_Pos[2] += dpos;
1507 temp_Pos[2] = tPos[2];
1510 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1511 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1512 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1520 double spacs[3], vec[3];
1521 double w, wx, wy, wz, len, fn, x, y, z, vol;
1522 double vol_density,sav;
1523 double *lower_corner, *upper_corner;
1532 for (i=0; i<3; i++) {
1533 len = upper_corner[i] - lower_corner[i];
1535 fn = len*vol_density + 1;
1536 npts[i] = (int)ceil(fn);
1537 spacs[i] = len/((double)(npts[i])-1.0);
1538 if (apolparm != VNULL) {
1540 if (apolparm->
grid[i] > spacs[i]) {
1541 Vnm_print(2,
"Vacc_totalSAV: Warning, your GRID value (%g) is larger than the recommended value (%g)!\n",
1542 apolparm->
grid[i], spacs[i]);
1544 spacs[i] = apolparm->
grid[i];
1550 for (x=lower_corner[0]; x<=upper_corner[0]; x=x+spacs[0]) {
1551 if ( VABS(x - lower_corner[0]) < VSMALL) {
1553 }
else if ( VABS(x - upper_corner[0]) < VSMALL) {
1559 for (y=lower_corner[1]; y<=upper_corner[1]; y=y+spacs[1]) {
1560 if ( VABS(y - lower_corner[1]) < VSMALL) {
1562 }
else if ( VABS(y - upper_corner[1]) < VSMALL) {
1568 for (z=lower_corner[2]; z<=upper_corner[2]; z=z+spacs[2]) {
1569 if ( VABS(z - lower_corner[2]) < VSMALL) {
1571 }
else if ( VABS(z - upper_corner[2]) < VSMALL) {
1586 w = spacs[0]*spacs[1]*spacs[2];
1593 Vclist *clist,
int iatom,
double *value) {
1599 int xmin, ymin, zmin;
1600 int xmax, ymax, zmax;
1602 double sigma6, sigma12;
1604 double spacs[3], vec[3];
1605 double w, wx, wy, wz, len, fn, x, y, z, vol;
1607 double vol_density, energy, rho,
srad;
1608 double psig, epsilon, watepsilon, sigma, watsigma, eni, chi;
1611 double *lower_corner, *upper_corner;
1613 Vatom *atom = VNULL;
1614 VASSERT(apolparm != VNULL);
1633 srad = apolparm->
srad;
1634 rho = apolparm->
bconc;
1639 sigma = psig + watsigma;
1640 epsilon = VSQRT((epsilon * watepsilon));
1643 sigma6 = VPOW(sigma,6);
1644 sigma12 = VPOW(sigma,12);
1647 xmin = pos[0] - pad;
1648 xmax = pos[0] + pad;
1649 ymin = pos[1] - pad;
1650 ymax = pos[1] + pad;
1651 zmin = pos[2] - pad;
1652 zmax = pos[2] + pad;
1654 for (i=0; i<3; i++) {
1655 len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1657 fn = len*vol_density + 1;
1658 npts[i] = (int)ceil(fn);
1661 if (apolparm->
grid[i] > spacs[i]) {
1662 Vnm_print(2,
"Vacc_totalSAV: Warning, your GRID value (%g) is larger than the recommended value (%g)!\n",
1663 apolparm->
grid[i], spacs[i]);
1665 spacs[i] = apolparm->
grid[i];
1669 for (x=xmin; x<=xmax; x=x+spacs[0]) {
1670 if ( VABS(x - xmin) < VSMALL) {
1672 }
else if ( VABS(x - xmax) < VSMALL) {
1678 for (y=ymin; y<=ymax; y=y+spacs[1]) {
1679 if ( VABS(y - ymin) < VSMALL) {
1681 }
else if ( VABS(y - ymax) < VSMALL) {
1687 for (z=zmin; z<=zmax; z=z+spacs[2]) {
1688 if ( VABS(z - zmin) < VSMALL) {
1690 }
else if ( VABS(z - zmax) < VSMALL) {
1701 if (VABS(chi) > VSMALL) {
1703 x2 = VSQR(vec[0]-pos[0]);
1704 y2 = VSQR(vec[1]-pos[1]);
1705 z2 = VSQR(vec[2]-pos[2]);
1706 r = VSQRT(x2+y2+z2);
1708 if (r <= 14 && r >= sigma) {
1709 eni = chi*rho*epsilon*(-2.0*sigma6/VPOW(r,6)+sigma12/VPOW(r,12));
1711 eni = -1.0*epsilon*chi*rho;
1725 w = spacs[0]*spacs[1]*spacs[2];
1739 double energy = 0.0;
1740 double tenergy = 0.0;
1741 double rho = apolparm->
bconc;
1745 if(apolparm->
setwat == 0){
1746 Vnm_print(2,
"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1750 if (VABS(rho) < VSMALL) {
1757 if(rc == 0)
return 0;
1822 VASSERT(apolparm != VNULL);
1826 if(apolparm->
setwat == 0){
1827 Vnm_print(2,
"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1839 srad = apolparm->
srad;
1840 rho = apolparm->
bconc;
1846 sigma = psig + watsigma;
1847 epsilon = VSQRT((epsilon * watepsilon));
1850 sigma6 = VPOW(sigma,6);
1851 sigma12 = VPOW(sigma,12);
1854 for (i=0; i<3; i++) {
1855 len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1857 fn = len*vol_density + 1;
1858 npts[i] = (int)ceil(fn);
1862 if (apolparm->
grid[i] > spacs[i]) {
1863 Vnm_print(2,
"Vacc_totalSAV: Warning, your GRID value (%g) is larger than the recommended value (%g)!\n",
1864 apolparm->
grid[i], spacs[i]);
1866 spacs[i] = apolparm->
grid[i];
1870 xmin = pos[0] - pad;
1871 xmax = pos[0] + pad;
1872 ymin = pos[1] - pad;
1873 ymax = pos[1] + pad;
1874 zmin = pos[2] - pad;
1875 zmax = pos[2] + pad;
1877 for (x=xmin; x<=xmax; x=x+spacs[0]) {
1878 if ( VABS(x - xmin) < VSMALL) {
1880 }
else if ( VABS(x - xmax) < VSMALL) {
1886 for (y=ymin; y<=ymax; y=y+spacs[1]) {
1887 if ( VABS(y - ymin) < VSMALL) {
1889 }
else if ( VABS(y - ymax) < VSMALL) {
1895 for (z=zmin; z<=zmax; z=z+spacs[2]) {
1896 if ( VABS(z - zmin) < VSMALL) {
1898 }
else if ( VABS(z - zmax) < VSMALL) {
1911 x2 = VSQR(vec[0]-pos[0]);
1912 y2 = VSQR(vec[1]-pos[1]);
1913 z2 = VSQR(vec[2]-pos[2]);
1914 r = VSQRT(x2+y2+z2);
1916 if (r <= 14 && r >= sigma){
1918 fo = 12.0*chi*rho*epsilon*(sigma6/VPOW(r,7)-sigma12/VPOW(r,13));
1920 fpt[0] = -1.0*(pos[0]-vec[0])*fo/r;
1921 fpt[1] = -1.0*(pos[1]-vec[1])*fo/r;
1922 fpt[2] = -1.0*(pos[2]-vec[2])*fo/r;
1925 for (si=0; si < 3; si++) fpt[si] = 0.0;
1928 for (si=0; si < 3; si++) fpt[si] = 0.0;
1932 force[i] += (w*fpt[i]);
1939 w = spacs[0]*spacs[1]*spacs[2];
1940 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.
#define MAX_SPHERE_PTS
Maximum number of points on a sphere.
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.