APBS 3.0.0
Loading...
Searching...
No Matches
vacc.c
Go to the documentation of this file.
1
57#include "vacc.h"
58
59VEMBED(rcsid="$Id$")
60
61#if !defined(VINLINE_VACC)
62
63VPUBLIC unsigned long int Vacc_memChk(Vacc *thee) {
64 if (thee == VNULL)
65 return 0;
66 return Vmem_bytes(thee->mem);
67}
68
69#endif /* if !defined(VINLINE_VACC) */
70
80VPRIVATE int ivdwAccExclus(
81 Vacc *thee,
82 double center[3],
83 double radius,
84 int atomID
85 ) {
86
87 int iatom;
88 double dist2,
89 *apos;
90 Vatom *atom;
91 VclistCell *cell;
92
93 VASSERT(thee != VNULL);
94
95 /* We can only test probes with radii less than the max specified */
96 if (radius > Vclist_maxRadius(thee->clist)) {
97 Vnm_print(2,
98 "Vacc_ivdwAcc: got radius (%g) bigger than max radius (%g)\n",
99 radius, Vclist_maxRadius(thee->clist));
100 VASSERT(0);
101 }
102
103 /* Get the relevant cell from the cell list */
104 cell = Vclist_getCell(thee->clist, center);
105
106 /* If we have no cell, then no atoms are nearby and we're definitely
107 * accessible */
108 if (cell == VNULL) {
109 return 1;
110 }
111
112 /* Otherwise, check for overlap with the atoms in the cell */
113 for (iatom=0; iatom<cell->natoms; iatom++) {
114 atom = cell->atoms[iatom];
115
116 // We don't actually need to test this if the atom IDs do match; don't compute this if we're comparing atom against itself.
117 if (atom->id == atomID) continue;
118
119 apos = atom->position;
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)){
123 return 0;
124 }
125 }
126
127 /* If we're still here, then the point is accessible */
128 return 1;
129
130}
131
132VPUBLIC Vacc* Vacc_ctor(Valist *alist,
133 Vclist *clist,
134 double surf_density /* Surface density */
135 ) {
136
137
138 Vacc *thee = VNULL;
139
140 /* Set up the structure */
141 thee = (Vacc*)Vmem_malloc(VNULL, 1, sizeof(Vacc) );
142 VASSERT( thee != VNULL);
143 VASSERT( Vacc_ctor2(thee, alist, clist, surf_density));
144 return thee;
145}
146
148VPRIVATE int Vacc_storeParms(Vacc *thee,
149 Valist *alist,
150 Vclist *clist,
151 double surf_density /* Surface density */
152 ) {
153
154 int nsphere,
155 iatom;
156 double maxrad = 0.0,
157 maxarea,
158 rad;
159 Vatom *atom;
160
161 if (alist == VNULL) {
162 Vnm_print(2, "Vacc_storeParms: Got NULL Valist!\n");
163 return 0;
164 } else thee->alist = alist;
165 if (clist == VNULL) {
166 Vnm_print(2, "Vacc_storeParms: Got NULL Vclist!\n");
167 return 0;
168 } else thee->clist = clist;
169 thee->surf_density = surf_density;
170
171 /* Loop through the atoms to determine the maximum radius */
172 maxrad = 0.0;
173 for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++) {
174 atom = Valist_getAtom(alist, iatom);
175 rad = Vatom_getRadius(atom);
176 if (rad > maxrad) maxrad = rad;
177 }
178 maxrad = maxrad + Vclist_maxRadius(thee->clist);
179
180 maxarea = 4.0*VPI*maxrad*maxrad;
181 nsphere = (int)ceil(maxarea*surf_density);
182
183 Vnm_print(0, "Vacc_storeParms: Surf. density = %g\n", surf_density);
184 Vnm_print(0, "Vacc_storeParms: Max area = %g\n", maxarea);
185 thee->refSphere = VaccSurf_refSphere(thee->mem, nsphere);
186 Vnm_print(0, "Vacc_storeParms: Using %d-point reference sphere\n",
187 thee->refSphere->npts);
188
189 return 1;
190}
191
193VPRIVATE int Vacc_allocate(Vacc *thee) {
194
195 int i,
196 natoms;
197
198 natoms = Valist_getNumberAtoms(thee->alist);
199
200 thee->atomFlags = (int*)Vmem_malloc(thee->mem, natoms, sizeof(int));
201 if (thee->atomFlags == VNULL) {
202 Vnm_print(2,
203 "Vacc_allocate: Failed to allocate %d (int)s for atomFlags!\n",
204 natoms);
205 return 0;
206 }
207 for (i=0; i<natoms; i++) (thee->atomFlags)[i] = 0;
208
209 return 1;
210}
211
212
213VPUBLIC int Vacc_ctor2(Vacc *thee,
214 Valist *alist,
215 Vclist *clist,
216 double surf_density
217 ) {
218
219 /* Check and store parameters */
220 if (!Vacc_storeParms(thee, alist, clist, surf_density)) {
221 Vnm_print(2, "Vacc_ctor2: parameter check failed!\n");
222 return 0;
223 }
224
225 /* Set up memory management object */
226 thee->mem = Vmem_ctor("APBS::VACC");
227 if (thee->mem == VNULL) {
228 Vnm_print(2, "Vacc_ctor2: memory object setup failed!\n");
229 return 0;
230 }
231
232 /* Setup and check probe */
233 thee->surf = VNULL;
234
235 /* Allocate space */
236 if (!Vacc_allocate(thee)) {
237 Vnm_print(2, "Vacc_ctor2: memory allocation failed!\n");
238 return 0;
239 }
240
241 return 1;
242}
243
244
245VPUBLIC void Vacc_dtor(Vacc **thee) {
246
247 if ((*thee) != VNULL) {
248 Vacc_dtor2(*thee);
249 Vmem_free(VNULL, 1, sizeof(Vacc), (void **)thee);
250 (*thee) = VNULL;
251 }
252
253}
254
255VPUBLIC void Vacc_dtor2(Vacc *thee) {
256
257 int i,
258 natoms;
259
260 natoms = Valist_getNumberAtoms(thee->alist);
261 Vmem_free(thee->mem, natoms, sizeof(int), (void **)&(thee->atomFlags));
262
263 if (thee->refSphere != VNULL) {
264 VaccSurf_dtor(&(thee->refSphere));
265 thee->refSphere = VNULL;
266 }
267 if (thee->surf != VNULL) {
268 for (i=0; i<natoms; i++) VaccSurf_dtor(&(thee->surf[i]));
269 Vmem_free(thee->mem, natoms, sizeof(VaccSurf *),
270 (void **)&(thee->surf));
271 thee->surf = VNULL;
272 }
273
274 Vmem_dtor(&(thee->mem));
275}
276
277VPUBLIC double Vacc_vdwAcc(Vacc *thee,
278 double center[3]
279 ) {
280
281 VclistCell *cell;
282 Vatom *atom;
283 int iatom;
284 double *apos,
285 dist2;
286
287 /* Get the relevant cell from the cell list */
288 cell = Vclist_getCell(thee->clist, center);
289
290 /* If we have no cell, then no atoms are nearby and we're definitely
291 * accessible */
292 if (cell == VNULL) return 1.0;
293
294 /* Otherwise, check for overlap with the atoms in the cell */
295 for (iatom=0; iatom<cell->natoms; iatom++) {
296 atom = cell->atoms[iatom];
297 apos = Vatom_getPosition(atom);
298 dist2 = VSQR(center[0]-apos[0]) + VSQR(center[1]-apos[1])
299 + VSQR(center[2]-apos[2]);
300 if (dist2 < VSQR(Vatom_getRadius(atom))) return 0.0;
301 }
302
303 /* If we're still here, then the point is accessible */
304 return 1.0;
305}
306
307VPUBLIC double Vacc_ivdwAcc(Vacc *thee,
308 double center[3],
309 double radius
310 ) {
311
312 return (double)ivdwAccExclus(thee, center, radius, -1);
313
314}
315
317 double center[VAPBS_DIM],
318 double win,
319 double infrad,
320 Vatom *atom,
321 double *grad
322 ) {
323
324 int i;
325 double dist,
326 *apos,
327 arad,
328 sm,
329 sm2,
330 w2i, /* inverse of win squared */
331 w3i, /* inverse of win cubed */
332 mygrad,
333 mychi = 1.0; /* Char. func. value for given atom */
334
335 VASSERT(thee != NULL);
336
337 /* Inverse squared window parameter */
338 w2i = 1.0/(win*win);
339 w3i = 1.0/(win*win*win);
340
341 /* The grad is zero by default */
342 for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
343
344 /* *** CALCULATE THE CHARACTERISTIC FUNCTION VALUE FOR THIS ATOM AND THE
345 * *** MAGNITUDE OF THE FORCE *** */
346 apos = Vatom_getPosition(atom);
347 /* Zero-radius atoms don't contribute */
348 if (Vatom_getRadius(atom) > 0.0) {
349 arad = Vatom_getRadius(atom) + infrad;
350 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
351 + VSQR(apos[2]-center[2]));
352 /* If we're inside an atom, the entire characteristic function
353 * will be zero and the grad will be zero, so we can stop */
354 if (dist < (arad - win)) return;
355 /* Likewise, if we're outside the smoothing window, the characteristic
356 * function is unity and the grad will be zero, so we can stop */
357 else if (dist > (arad + win)) return;
358 /* Account for floating point error at the border
359 * NAB: COULDN'T THESE TESTS BE COMBINED AS BELOW
360 * (Vacc_splineAccAtom)? */
361 else if ((VABS(dist - (arad - win)) < VSMALL) ||
362 (VABS(dist - (arad + win)) < VSMALL)) return;
363 /* If we're inside the smoothing window */
364 else {
365 sm = dist - arad + win;
366 sm2 = VSQR(sm);
367 mychi = 0.75*sm2*w2i -0.25*sm*sm2*w3i;
368 mygrad = 1.5*sm*w2i - 0.75*sm2*w3i;
369 }
370 /* Now assemble the grad vector */
371 VASSERT(mychi > 0.0);
372 for (i=0; i<VAPBS_DIM; i++)
373 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
374 }
375}
376
378 double center[VAPBS_DIM],
379 double win,
380 double infrad,
381 Vatom *atom,
382 double *grad
383 ) {
384
385 int i;
386 double dist,
387 *apos,
388 arad,
389 sm,
390 sm2,
391 w2i, /* Inverse of win squared */
392 w3i, /* Inverse of win cubed */
393 mygrad,
394 mychi = 1.0; /* Char. func. value for given atom */
395
396 VASSERT(thee != NULL);
397
398 /* Inverse squared window parameter */
399 w2i = 1.0/(win*win);
400 w3i = 1.0/(win*win*win);
401
402 /* The grad is zero by default */
403 for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
404
405 /* *** CALCULATE THE CHARACTERISTIC FUNCTION VALUE FOR THIS ATOM AND THE
406 * *** MAGNITUDE OF THE FORCE *** */
407 apos = Vatom_getPosition(atom);
408 /* Zero-radius atoms don't contribute */
409 if (Vatom_getRadius(atom) > 0.0) {
410 arad = Vatom_getRadius(atom) + infrad;
411 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
412 + VSQR(apos[2]-center[2]));
413 /* If we're inside an atom, the entire characteristic function
414 * will be zero and the grad will be zero, so we can stop */
415 if (dist < (arad - win)) return;
416 /* Likewise, if we're outside the smoothing window, the characteristic
417 * function is unity and the grad will be zero, so we can stop */
418 else if (dist > (arad + win)) return;
419 /* Account for floating point error at the border
420 * NAB: COULDN'T THESE TESTS BE COMBINED AS BELOW
421 * (Vacc_splineAccAtom)? */
422 else if ((VABS(dist - (arad - win)) < VSMALL) ||
423 (VABS(dist - (arad + win)) < VSMALL)) return;
424 /* If we're inside the smoothing window */
425 else {
426 sm = dist - arad + win;
427 sm2 = VSQR(sm);
428 mychi = 0.75*sm2*w2i -0.25*sm*sm2*w3i;
429 mygrad = 1.5*sm*w2i - 0.75*sm2*w3i;
430 }
431 /* Now assemble the grad vector */
432 VASSERT(mychi > 0.0);
433 for (i=0; i<VAPBS_DIM; i++)
434 grad[i] = -(mygrad)*((center[i] - apos[i])/dist);
435 }
436}
437
438VPUBLIC double Vacc_splineAccAtom(Vacc *thee,
439 double center[VAPBS_DIM],
440 double win,
441 double infrad,
442 Vatom *atom
443 ) {
444
445 double dist,
446 *apos,
447 arad,
448 sm,
449 sm2,
450 w2i, /* Inverse of win squared */
451 w3i, /* Inverse of win cubed */
452 value,
453 stot,
454 sctot;
455
456 VASSERT(thee != NULL);
457
458 /* Inverse squared window parameter */
459 w2i = 1.0/(win*win);
460 w3i = 1.0/(win*win*win);
461
462 apos = Vatom_getPosition(atom);
463 /* Zero-radius atoms don't contribute */
464 if (Vatom_getRadius(atom) > 0.0) {
465 arad = Vatom_getRadius(atom) + infrad;
466 stot = arad + 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]));
470 /* If we're inside an atom, the entire characteristic function
471 * will be zero */
472 if ((dist < sctot) || (VABS(dist - sctot) < VSMALL)){
473 value = 0.0;
474 /* We're outside the smoothing window */
475 } else if ((dist > stot) || (VABS(dist - stot) < VSMALL)) {
476 value = 1.0;
477 /* We're inside the smoothing window */
478 } else {
479 sm = dist - arad + win;
480 sm2 = VSQR(sm);
481 value = 0.75*sm2*w2i - 0.25*sm*sm2*w3i;
482 }
483 } else value = 1.0;
484
485 return value;
486}
487
493VPRIVATE double splineAcc(
494 Vacc *thee,
495 double center[VAPBS_DIM],
497 double win,
498 double infrad,
499 VclistCell *cell
500 ) {
501
502 int atomID, iatom;
503 Vatom *atom;
504 double value = 1.0;
505
506 VASSERT(thee != NULL);
507
508 /* Now loop through the atoms assembling the characteristic function */
509 for (iatom=0; iatom<cell->natoms; iatom++) {
510
511 atom = cell->atoms[iatom];
512 atomID = atom->id;
513
514 /* Check to see if we've counted this atom already */
515 if ( !(thee->atomFlags[atomID]) ) {
516
517 thee->atomFlags[atomID] = 1;
518 value *= Vacc_splineAccAtom(thee, center, win, infrad, atom);
519
520 if (value < VSMALL) return value;
521 }
522 }
523
524 return value;
525}
526
527
528VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win,
529 double infrad) {
530
531 VclistCell *cell;
532 Vatom *atom;
533 int iatom, atomID;
534
535
536 VASSERT(thee != NULL);
537
538 if (Vclist_maxRadius(thee->clist) < (win + infrad)) {
539 Vnm_print(2, "Vacc_splineAcc: Vclist has max_radius=%g;\n",
540 Vclist_maxRadius(thee->clist));
541 Vnm_print(2, "Vacc_splineAcc: Insufficient for win=%g, infrad=%g\n",
542 win, infrad);
543 VASSERT(0);
544 }
545
546 /* Get a cell or VNULL; in the latter case return 1.0 */
547 cell = Vclist_getCell(thee->clist, center);
548 if (cell == VNULL) return 1.0;
549
550 /* First, reset the list of atom flags
551 * NAB: THIS SEEMS VERY INEFFICIENT */
552 for (iatom=0; iatom<cell->natoms; iatom++) {
553 atom = cell->atoms[iatom];
554 atomID = atom->id;
555 thee->atomFlags[atomID] = 0;
556 }
557
558 return splineAcc(thee, center, win, infrad, cell);
559}
560
561VPUBLIC void Vacc_splineAccGrad(Vacc *thee, double center[VAPBS_DIM],
562 double win, double infrad, double *grad) {
563
564 int iatom, i, atomID;
565 double acc = 1.0;
566 double tgrad[VAPBS_DIM];
567 VclistCell *cell;
568 Vatom *atom = VNULL;
569
570 VASSERT(thee != NULL);
571
572 if (Vclist_maxRadius(thee->clist) < (win + infrad)) {
573 Vnm_print(2, "Vacc_splineAccGrad: Vclist max_radius=%g;\n",
574 Vclist_maxRadius(thee->clist));
575 Vnm_print(2, "Vacc_splineAccGrad: Insufficient for win=%g, infrad=%g\n",
576 win, infrad);
577 VASSERT(0);
578 }
579
580 /* Reset the gradient */
581 for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
582
583 /* Get the cell; check for nullity */
584 cell = Vclist_getCell(thee->clist, center);
585 if (cell == VNULL) return;
586
587 /* Reset the list of atom flags */
588 for (iatom=0; iatom<cell->natoms; iatom++) {
589 atom = cell->atoms[iatom];
590 atomID = atom->id;
591 thee->atomFlags[atomID] = 0;
592 }
593
594 /* Get the local accessibility */
595 acc = splineAcc(thee, center, win, infrad, cell);
596
597 /* Accumulate the gradient of all local atoms */
598 if (acc > VSMALL) {
599 for (iatom=0; iatom<cell->natoms; iatom++) {
600 atom = cell->atoms[iatom];
601 Vacc_splineAccGradAtomNorm(thee, center, win, infrad, atom, tgrad);
602 }
603 for (i=0; i<VAPBS_DIM; i++) grad[i] += tgrad[i];
604 }
605 for (i=0; i<VAPBS_DIM; i++) grad[i] *= -acc;
606}
607
608VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM],
609 double radius) {
610
611 double rc;
612
613 /* ******* CHECK IF OUTSIDE ATOM+PROBE RADIUS SURFACE ***** */
614 if (Vacc_ivdwAcc(thee, center, radius) == 1.0) {
615
616 /* Vnm_print(2, "DEBUG: ivdwAcc = 1.0\n"); */
617 rc = 1.0;
618
619 /* ******* CHECK IF INSIDE ATOM RADIUS SURFACE ***** */
620 } else if (Vacc_vdwAcc(thee, center) == 0.0) {
621
622 /* Vnm_print(2, "DEBUG: vdwAcc = 0.0\n"); */
623 rc = 0.0;
624
625 /* ******* CHECK IF OUTSIDE MOLECULAR SURFACE ***** */
626 } else {
627
628 /* Vnm_print(2, "DEBUG: calling fastMolAcc...\n"); */
629 rc = Vacc_fastMolAcc(thee, center, radius);
630
631 }
632
633 return rc;
634
635}
636
637VPUBLIC double Vacc_fastMolAcc(Vacc *thee, double center[VAPBS_DIM],
638 double radius) {
639
640 Vatom *atom;
641 VaccSurf *surf;
642 VclistCell *cell;
643 int ipt, iatom, atomID;
644 double dist2, rad2;
645
646 rad2 = radius*radius;
647
648 /* Check to see if the SAS has been defined */
649 if (thee->surf == VNULL) Vacc_SASA(thee, radius);
650
651 /* Get the cell associated with this point */
652 cell = Vclist_getCell(thee->clist, center);
653 if (cell == VNULL) {
654 Vnm_print(2, "Vacc_fastMolAcc: unexpected VNULL VclistCell!\n");
655 return 1.0;
656 }
657
658 /* Loop through all the atoms in the cell */
659 for (iatom=0; iatom<cell->natoms; iatom++) {
660 atom = cell->atoms[iatom];
661 atomID = Vatom_getAtomID(atom);
662 surf = thee->surf[atomID];
663 /* Loop through all SAS points associated with this atom */
664 for (ipt=0; ipt<surf->npts; ipt++) {
665 /* See if we're within a probe radius of the point */
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;
670 }
671 }
672
673 /* If all else failed, we are not inside the molecular surface */
674 return 0.0;
675}
676
677
678#if defined(HAVE_MC_H)
679VPUBLIC void Vacc_writeGMV(Vacc *thee, double radius, int meth, Gem *gm,
680 char *iodev, char *iofmt, char *iohost, char *iofile) {
681
682 double *accVals[MAXV], coord[3];
683 Vio *sock;
684 int ivert, icoord;
685
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);
692 if (meth == 0) {
693 accVals[0][ivert] = Vacc_molAcc(thee, coord, radius);
694 accVals[1][ivert] = Vacc_molAcc(thee, coord, radius);
695 } else if (meth == 1) {
696 accVals[0][ivert] = Vacc_ivdwAcc(thee, coord, radius);
697 accVals[1][ivert] = Vacc_ivdwAcc(thee, coord, radius);
698 } else if (meth == 2) {
699 accVals[0][ivert] = Vacc_vdwAcc(thee, coord);
700 accVals[1][ivert] = Vacc_vdwAcc(thee, coord);
701 } else VASSERT(0);
702 }
703 sock = Vio_ctor(iodev, iofmt, iohost, iofile, "w");
704 Gem_writeGMV(gm, sock, 1, accVals);
705 Vio_dtor(&sock);
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]));
710}
711#endif /* defined(HAVE_MC_H) */
712
713VPUBLIC double Vacc_SASA(Vacc *thee,
714 double radius
715 ) {
716
717 int i,
718 natom;
719 double area;
720 //*apos; // gcc says unused
721 Vatom *atom;
722 VaccSurf *asurf;
723
724 time_t ts; // PCE: temp
725 ts = clock();
726
727 //unsigned long long mbeg; // gcc says unused
728
729 natom = Valist_getNumberAtoms(thee->alist);
730
731 /* Check to see if we need to build the surface */
732 if (thee->surf == VNULL) {
733 thee->surf = Vmem_malloc(thee->mem, natom, sizeof(VaccSurf *));
734
735#if defined(DEBUG_MAC_OSX_OCL) || defined(DEBUG_MAC_OSX_STANDARD)
736#include "mach_chud.h"
737 machm_(&mbeg);
738#pragma omp parallel for private(i,atom)
739#endif
740 for (i=0; i<natom; i++) {
741 atom = Valist_getAtom(thee->alist, i);
742 /* NOTE: RIGHT NOW WE DO THIS FOR THE ENTIRE MOLECULE WHICH IS
743 * INCREDIBLY INEFFICIENT, PARTICULARLY DURING FOCUSING!!! */
744 thee->surf[i] = Vacc_atomSurf(thee, atom, thee->refSphere,
745 radius);
746 }
747 }
748
749 /* Calculate the area */
750 area = 0.0;
751 for (i=0; i<natom; i++) {
752 atom = Valist_getAtom(thee->alist, i);
753 asurf = thee->surf[i];
754 /* See if this surface needs to be rebuilt */
755 if (asurf->probe_radius != radius) {
756 Vnm_print(2, "Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
757 asurf->probe_radius, radius);
758 VaccSurf_dtor2(asurf);
759 thee->surf[i] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
760 asurf = thee->surf[i];
761 }
762 area += (asurf->area);
763 }
764
765#if defined(DEBUG_MAC_OSX_OCL) || defined(DEBUG_MAC_OSX_STANDARD)
766 mets_(&mbeg, "Vacc_SASA - Parallel");
767#endif
768
769 Vnm_print(0, "Vacc_SASA: Time elapsed: %f\n", ((double)clock() - ts) / CLOCKS_PER_SEC);
770 return area;
771
772}
773
774VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius) {
775
776 return Vacc_SASA(thee, radius);
777
778}
779
780VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom) {
781
782 VaccSurf *asurf;
783 int id;
784
785 if (thee->surf == VNULL) Vacc_SASA(thee, radius);
786
787 id = Vatom_getAtomID(atom);
788 asurf = thee->surf[id];
789
790 /* See if this surface needs to be rebuilt */
791 if (asurf->probe_radius != radius) {
792 Vnm_print(2, "Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
793 asurf->probe_radius, radius);
794 VaccSurf_dtor2(asurf);
795 thee->surf[id] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
796 asurf = thee->surf[id];
797 }
798
799 return asurf->area;
800
801}
802
803VPUBLIC VaccSurf* VaccSurf_ctor(Vmem *mem, double probe_radius, int nsphere) {
804 VaccSurf *thee;
805
806 //thee = Vmem_malloc(mem, 1, sizeof(Vacc) );
807 thee = (VaccSurf*)calloc(1,sizeof(Vacc));
808 VASSERT( VaccSurf_ctor2(thee, mem, probe_radius, nsphere) );
809
810 return thee;
811}
812
813VPUBLIC int VaccSurf_ctor2(VaccSurf *thee, Vmem *mem, double probe_radius,
814 int nsphere) {
815
816 if (thee == VNULL)
817 return 0;
818
819 thee->mem = mem;
820 thee->npts = nsphere;
821 thee->probe_radius = probe_radius;
822 thee->area = 0.0;
823
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));
829 } else {
830 thee->xpts = VNULL;
831 thee->ypts = VNULL;
832 thee->zpts = VNULL;
833 thee->bpts = VNULL;
834 }
835
836 return 1;
837}
838
839VPUBLIC void VaccSurf_dtor(VaccSurf **thee) {
840
841 Vmem *mem;
842
843 if ((*thee) != VNULL) {
844 mem = (*thee)->mem;
845 VaccSurf_dtor2(*thee);
846 //Vmem_free(mem, 1, sizeof(VaccSurf), (void **)thee);
847 free(*thee);
848 (*thee) = VNULL;
849 }
850
851}
852
853VPUBLIC void VaccSurf_dtor2(VaccSurf *thee) {
854
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));
864 }
865
866}
867
868VPUBLIC VaccSurf* Vacc_atomSurf(Vacc *thee, Vatom *atom,
869 VaccSurf *ref, double prad) {
870
871 VaccSurf *surf;
872 size_t i, j, npts;
873 int atomID;
874 double arad, rad, pos[3], *apos;
875
876 /* Get atom information */
877 arad = Vatom_getRadius(atom);
878 apos = Vatom_getPosition(atom);
879 atomID = Vatom_getAtomID(atom);
880
881 if (arad < VSMALL) {
882 return VaccSurf_ctor(thee->mem, prad, 0);
883 }
884
885 rad = arad + prad;
886
887 /* Determine which points will contribute */
888 npts = 0;
889 for (i=0; i<ref->npts; i++) {
890 /* Reset point flag: zero-radius atoms do not contribute */
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];
894 if (ivdwAccExclus(thee, pos, prad, atomID)) {
895 npts++;
896 ref->bpts[i] = (ref->bpts[i] << 1) + 1;
897 } else {
898 ref->bpts[i] <<= 1;
899 }
900 }
901
902 /* Allocate space for the points */
903 surf = VaccSurf_ctor(thee->mem, prad, npts);
904
905 /* Assign the points */
906 j = 0;
907 for (i=0; i<ref->npts; i++) {
908 char flag = ref->bpts[i] & 1;
909 ref->bpts[i] >>= 1;
910 if (flag) {
911 surf->bpts[j] = 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];
915 j++;
916 }
917 }
918
919 /* Assign the area */
920 surf->area = 4.0*VPI*rad*rad*((double)(surf->npts))/((double)(ref->npts));
921
922 return surf;
923
924}
925
926VPUBLIC VaccSurf* VaccSurf_refSphere(Vmem *mem, int npts) {
927
928 VaccSurf *surf;
929 int nactual, i, itheta, ntheta, iphi, nphimax, nphi;
930 double frac;
931 double sintheta, costheta, theta, dtheta;
932 double sinphi, cosphi, phi, dphi;
933
934 /* Setup "constants" */
935 frac = ((double)(npts))/4.0;
936 ntheta = VRINT(VSQRT(Vunit_pi*frac));
937 dtheta = Vunit_pi/((double)(ntheta));
938 nphimax = 2*ntheta;
939
940 /* Count the actual number of points to be used */
941 nactual = 0;
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);
947 nactual += nphi;
948 }
949
950 /* Allocate space for the points */
951 surf = VaccSurf_ctor(mem, 1.0, nactual);
952
953 /* Clear out the boolean array */
954 for (i=0; i<nactual; i++) surf->bpts[i] = 1;
955
956 /* Assign the points */
957 nactual = 0;
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);
963 if (nphi != 0) {
964 dphi = 2*Vunit_pi/((double)(nphi));
965 for (iphi=0; iphi<nphi; iphi++) {
966 phi = dphi*((double)(iphi));
967 sinphi = VSIN(phi);
968 cosphi = VCOS(phi);
969 surf->xpts[nactual] = cosphi * sintheta;
970 surf->ypts[nactual] = sinphi * sintheta;
971 surf->zpts[nactual] = costheta;
972 nactual++;
973 }
974 }
975 }
976
977 surf->npts = nactual;
978
979 return surf;
980}
981
982VPUBLIC VaccSurf* Vacc_atomSASPoints(Vacc *thee, double radius,
983 Vatom *atom) {
984
985 VaccSurf *asurf = VNULL;
986 int id;
987
988 if (thee->surf == VNULL) Vacc_SASA(thee, radius);
989 id = Vatom_getAtomID(atom);
990
991 asurf = thee->surf[id];
992
993 /* See if this surface needs to be rebuilt */
994 if (asurf->probe_radius != radius) {
995 Vnm_print(2, "Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
996 asurf->probe_radius, radius);
997 VaccSurf_dtor2(asurf);
998 thee->surf[id] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
999 asurf = thee->surf[id];
1000 }
1001
1002 return asurf;
1003
1004}
1005
1006VPUBLIC void Vacc_splineAccGradAtomNorm4(Vacc *thee, double center[VAPBS_DIM],
1007 double win, double infrad, Vatom *atom, double *grad) {
1008
1009 int i;
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;
1015 double mychi = 1.0; /* Char. func. value for given atom */
1016
1017 VASSERT(thee != NULL);
1018
1019 /* The grad is zero by default */
1020 for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
1021
1022 /* *** CALCULATE THE CHARACTERISTIC FUNCTION VALUE FOR THIS ATOM AND THE
1023 * *** MAGNITUDE OF THE FORCE *** */
1024 apos = Vatom_getPosition(atom);
1025 /* Zero-radius atoms don't contribute */
1026 if (Vatom_getRadius(atom) > 0.0) {
1027
1028 arad = Vatom_getRadius(atom);
1029 arad = arad + infrad;
1030 b = arad - win;
1031 e = arad + win;
1032
1033 e2 = e * e;
1034 e3 = e2 * e;
1035 e4 = e3 * e;
1036 e5 = e4 * e;
1037 e6 = e5 * e;
1038 e7 = e6 * e;
1039 b2 = b * b;
1040 b3 = b2 * b;
1041 b4 = b3 * b;
1042 b5 = b4 * b;
1043 b6 = b5 * b;
1044 b7 = b6 * b;
1045
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;
1055 c7 = -20.0/denom;
1056
1057 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1058 + VSQR(apos[2]-center[2]));
1059
1060 /* If we're inside an atom, the entire characteristic function
1061 * will be zero and the grad will be zero, so we can stop */
1062 if (dist < (arad - win)) return;
1063 /* Likewise, if we're outside the smoothing window, the characteristic
1064 * function is unity and the grad will be zero, so we can stop */
1065 else if (dist > (arad + win)) return;
1066 /* Account for floating point error at the border
1067 * NAB: COULDN'T THESE TESTS BE COMBINED AS BELOW
1068 * (Vacc_splineAccAtom)? */
1069 else if ((VABS(dist - (arad - win)) < VSMALL) ||
1070 (VABS(dist - (arad + win)) < VSMALL)) return;
1071 /* If we're inside the smoothing window */
1072 else {
1073 sm = dist;
1074 sm2 = sm * sm;
1075 sm3 = sm2 * sm;
1076 sm4 = sm3 * sm;
1077 sm5 = sm4 * sm;
1078 sm6 = sm5 * sm;
1079 sm7 = sm6 * sm;
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;
1084 if (mychi <= 0.0) {
1085 /* Avoid numerical round off errors */
1086 return;
1087 } else if (mychi > 1.0) {
1088 /* Avoid numerical round off errors */
1089 mychi = 1.0;
1090 }
1091 }
1092 /* Now assemble the grad vector */
1093 VASSERT(mychi > 0.0);
1094 for (i=0; i<VAPBS_DIM; i++)
1095 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1096 }
1097}
1098
1099VPUBLIC void Vacc_splineAccGradAtomNorm3(Vacc *thee, double center[VAPBS_DIM],
1100 double win, double infrad, Vatom *atom, double *grad) {
1101
1102 int i;
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;
1108 double mychi = 1.0; /* Char. func. value for given atom */
1109
1110 VASSERT(thee != NULL);
1111
1112 /* The grad is zero by default */
1113 for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
1114
1115 /* *** CALCULATE THE CHARACTERISTIC FUNCTION VALUE FOR THIS ATOM AND THE
1116 * *** MAGNITUDE OF THE FORCE *** */
1117 apos = Vatom_getPosition(atom);
1118 /* Zero-radius atoms don't contribute */
1119 if (Vatom_getRadius(atom) > 0.0) {
1120
1121 arad = Vatom_getRadius(atom);
1122 arad = arad + infrad;
1123 b = arad - win;
1124 e = arad + win;
1125
1126 e2 = e * e;
1127 e3 = e2 * e;
1128 e4 = e3 * e;
1129 e5 = e4 * e;
1130 b2 = b * b;
1131 b3 = b2 * b;
1132 b4 = b3 * b;
1133 b5 = b4 * b;
1134
1135 denom = pow((e - b), 5.0);
1136 c0 = -10.0*e2*b3 + 5.0*e*b4 - b5;
1137 c1 = 30.0*e2*b2;
1138 c2 = -30.0*(e2*b + e*b2);
1139 c3 = 10.0*(e2 + 4.0*e*b + b2);
1140 c4 = -15.0*(e + b);
1141 c5 = 6;
1142 c0 = c0/denom;
1143 c1 = c1/denom;
1144 c2 = c2/denom;
1145 c3 = c3/denom;
1146 c4 = c4/denom;
1147 c5 = c5/denom;
1148
1149 dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1150 + VSQR(apos[2]-center[2]));
1151
1152 /* If we're inside an atom, the entire characteristic function
1153 * will be zero and the grad will be zero, so we can stop */
1154 if (dist < (arad - win)) return;
1155 /* Likewise, if we're outside the smoothing window, the characteristic
1156 * function is unity and the grad will be zero, so we can stop */
1157 else if (dist > (arad + win)) return;
1158 /* Account for floating point error at the border
1159 * NAB: COULDN'T THESE TESTS BE COMBINED AS BELOW
1160 * (Vacc_splineAccAtom)? */
1161 else if ((VABS(dist - (arad - win)) < VSMALL) ||
1162 (VABS(dist - (arad + win)) < VSMALL)) return;
1163 /* If we're inside the smoothing window */
1164 else {
1165 sm = dist;
1166 sm2 = sm * sm;
1167 sm3 = sm2 * sm;
1168 sm4 = sm3 * sm;
1169 sm5 = sm4 * sm;
1170 mychi = c0 + c1*sm + c2*sm2 + c3*sm3
1171 + c4*sm4 + c5*sm5;
1172 mygrad = c1 + 2.0*c2*sm + 3.0*c3*sm2 + 4.0*c4*sm3
1173 + 5.0*c5*sm4;
1174 if (mychi <= 0.0) {
1175 /* Avoid numerical round off errors */
1176 return;
1177 } else if (mychi > 1.0) {
1178 /* Avoid numerical round off errors */
1179 mychi = 1.0;
1180 }
1181 }
1182 /* Now assemble the grad vector */
1183 VASSERT(mychi > 0.0);
1184 for (i=0; i<VAPBS_DIM; i++)
1185 grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1186 }
1187}
1188
1189/* ///////////////////////////////////////////////////////////////////////////
1190 // Routine: Vacc_atomdSAV
1191 //
1192 // Purpose: Calculates the vector valued atomic derivative of volume
1193 //
1194 // Args: radius The radius of the solvent probe in Angstroms
1195 // iatom Index of the atom in thee->alist
1196 //
1197 // Author: Jason Wagoner
1198 // Nathan Baker (original FORTRAN routine from UHBD by Brock Luty)
1200VPUBLIC void Vacc_atomdSAV(Vacc *thee,
1201 double srad,
1202 Vatom *atom,
1203 double *dSA
1204 ) {
1205
1206 int ipt, iatom;
1207
1208 double area;
1209 double *tPos, tRad, vec[3];
1210 double dx,dy,dz;
1211 VaccSurf *ref;
1212 dx = 0.0;
1213 dy = 0.0;
1214 dz = 0.0;
1215 /* Get the atom information */
1216 ref = thee->refSphere;
1217 iatom = Vatom_getAtomID(atom);
1218
1219 dSA[0] = 0.0;
1220 dSA[1] = 0.0;
1221 dSA[2] = 0.0;
1222
1223 tPos = Vatom_getPosition(atom);
1224 tRad = Vatom_getRadius(atom);
1225
1226 if(tRad == 0.0) return;
1227
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];
1233 if (ivdwAccExclus(thee, vec, srad, iatom)) {
1234 dx = dx+vec[0]-tPos[0];
1235 dy = dy+vec[1]-tPos[1];
1236 dz = dz+vec[2]-tPos[2];
1237 }
1238 }
1239
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);
1244 }
1245
1246}
1247
1248/* Note: This is purely test code to make certain that the dSASA code is
1249 behaving properly. This function should NEVER be called by anyone
1250 other than an APBS developer at Wash U.
1251*/
1252VPRIVATE double Vacc_SASAPos(Vacc *thee, double radius) {
1253
1254 int i, natom;
1255 double area;
1256 Vatom *atom;
1257 VaccSurf *asurf;
1258
1259 natom = Valist_getNumberAtoms(thee->alist);
1260
1261 /* Calculate the area */
1262 area = 0.0;
1263 for (i=0; i<natom; i++) {
1264 atom = Valist_getAtom(thee->alist, i);
1265 asurf = thee->surf[i];
1266
1267 VaccSurf_dtor2(asurf);
1268 thee->surf[i] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
1269 asurf = thee->surf[i];
1270 area += (asurf->area);
1271 }
1272
1273 return area;
1274
1275}
1276
1277VPRIVATE double Vacc_atomSASAPos(Vacc *thee,
1278 double radius,
1279 Vatom *atom, /* The atom being manipulated */
1280 int mode
1281 ) {
1282
1283 VaccSurf *asurf;
1284 int id;
1285 static int warned = 0;
1286
1287 if ((thee->surf == VNULL) || (mode == 1)){
1288 if(!warned){
1289 Vnm_print(2, "WARNING: Recalculating entire surface!!!!\n");
1290 warned = 1;
1291 }
1292 Vacc_SASAPos(thee, radius); // reinitialize before we can do anything about doing a calculation on a repositioned atom
1293 }
1294
1295 id = Vatom_getAtomID(atom);
1296 asurf = thee->surf[id];
1297
1298 VaccSurf_dtor(&asurf);
1299 thee->surf[id] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
1300 asurf = thee->surf[id];
1301
1302 //printf("%s: Time elapsed: %f\n", __func__, ((double)clock() - ts) / CLOCKS_PER_SEC);
1303
1304 return asurf->area;
1305}
1306
1307/* ///////////////////////////////////////////////////////////////////////////
1308 // Routine: Vacc_atomdSASA
1309 //
1310 // Purpose: Calculates the derivative of surface area with respect to atomic
1311 // displacement using finite difference methods.
1312 //
1313 // Args: radius The radius of the solvent probe in Angstroms
1314 // iatom Index of the atom in thee->alist
1315 //
1316 // Author: Jason Wagoner
1317 // David Gohara
1318 // Nathan Baker (original FORTRAN routine from UHBD by Brock Luty)
1320VPUBLIC void Vacc_atomdSASA(Vacc *thee,
1321 double dpos,
1322 double srad,
1323 Vatom *atom,
1324 double *dSA
1325 ) {
1326
1327 double *temp_Pos,
1328 tPos[3],
1329 axb1,
1330 axt1,
1331 ayb1,
1332 ayt1,
1333 azb1,
1334 azt1;
1335 VaccSurf *ref;
1336
1337 //printf("%s: entering\n", __func__);
1338 time_t ts;
1339 ts = clock();
1340
1341 /* Get the atom information */
1342 ref = thee->refSphere;
1343 temp_Pos = Vatom_getPosition(atom); // Get a pointer to the position object. You actually manipulate the atom doing this...
1344
1345 tPos[0] = temp_Pos[0];
1346 tPos[1] = temp_Pos[1];
1347 tPos[2] = temp_Pos[2];
1348
1349 /* Shift by pos -/+ on x */
1350 temp_Pos[0] -= dpos;
1351 axb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1352 temp_Pos[0] = tPos[0];
1353
1354 temp_Pos[0] += dpos;
1355 axt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1356 temp_Pos[0] = tPos[0];
1357
1358 /* Shift by pos -/+ on y */
1359 temp_Pos[1] -= dpos;
1360 ayb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1361 temp_Pos[1] = tPos[1];
1362
1363 temp_Pos[1] += dpos;
1364 ayt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1365 temp_Pos[1] = tPos[1];
1366
1367 /* Shift by pos -/+ on z */
1368 temp_Pos[2] -= dpos;
1369 azb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1370 temp_Pos[2] = tPos[2];
1371
1372 temp_Pos[2] += dpos;
1373 azt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1374 temp_Pos[2] = tPos[2];
1375
1376 /* Reset the atom SASA to zero displacement */
1377 Vacc_atomSASAPos(thee, srad, atom,0);
1378
1379 /* Calculate the final value */
1380 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1381 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1382 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1383}
1384
1385/* Note: This is purely test code to make certain that the dSASA code is
1386 behaving properly. This function should NEVER be called by anyone
1387 other than an APBS developer at Wash U.
1388*/
1389VPUBLIC void Vacc_totalAtomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA) {
1390
1391 int iatom;
1392 double *temp_Pos, tRad;
1393 double tPos[3];
1394 double axb1,axt1,ayb1,ayt1,azb1,azt1;
1395 VaccSurf *ref;
1396
1397 /* Get the atom information */
1398 ref = thee->refSphere;
1399 temp_Pos = Vatom_getPosition(atom);
1400 tRad = Vatom_getRadius(atom);
1401 iatom = Vatom_getAtomID(atom);
1402
1403 dSA[0] = 0.0;
1404 dSA[1] = 0.0;
1405 dSA[2] = 0.0;
1406
1407 tPos[0] = temp_Pos[0];
1408 tPos[1] = temp_Pos[1];
1409 tPos[2] = temp_Pos[2];
1410
1411 /* Shift by pos -/+ on x */
1412 temp_Pos[0] -= dpos;
1413 axb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1414 temp_Pos[0] = tPos[0];
1415
1416 temp_Pos[0] += dpos;
1417 axt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1418 temp_Pos[0] = tPos[0];
1419
1420 /* Shift by pos -/+ on y */
1421 temp_Pos[1] -= dpos;
1422 ayb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1423 temp_Pos[1] = tPos[1];
1424
1425 temp_Pos[1] += dpos;
1426 ayt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1427 temp_Pos[1] = tPos[1];
1428
1429 /* Shift by pos -/+ on z */
1430 temp_Pos[2] -= dpos;
1431 azb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1432 temp_Pos[2] = tPos[2];
1433
1434 temp_Pos[2] += dpos;
1435 azt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1436 temp_Pos[2] = tPos[2];
1437
1438 /* Calculate the final value */
1439 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1440 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1441 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1442}
1443
1444/* Note: This is purely test code to make certain that the dSASA code is
1445 behaving properly. This function should NEVER be called by anyone
1446 other than an APBS developer at Wash U.
1447*/
1448VPUBLIC void Vacc_totalAtomdSAV(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA, Vclist *clist) {
1449
1450 int iatom;
1451 double *temp_Pos, tRad;
1452 double tPos[3];
1453 double axb1,axt1,ayb1,ayt1,azb1,azt1;
1454 VaccSurf *ref;
1455
1456 /* Get the atom information */
1457 ref = thee->refSphere;
1458 temp_Pos = Vatom_getPosition(atom);
1459 tRad = Vatom_getRadius(atom);
1460 iatom = Vatom_getAtomID(atom);
1461
1462 dSA[0] = 0.0;
1463 dSA[1] = 0.0;
1464 dSA[2] = 0.0;
1465
1466 tPos[0] = temp_Pos[0];
1467 tPos[1] = temp_Pos[1];
1468 tPos[2] = temp_Pos[2];
1469
1470 /* Shift by pos -/+ on x */
1471 temp_Pos[0] -= dpos;
1472 axb1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1473 temp_Pos[0] = tPos[0];
1474
1475 temp_Pos[0] += dpos;
1476 axt1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1477 temp_Pos[0] = tPos[0];
1478
1479 /* Shift by pos -/+ on y */
1480 temp_Pos[1] -= dpos;
1481 ayb1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1482 temp_Pos[1] = tPos[1];
1483
1484 temp_Pos[1] += dpos;
1485 ayt1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1486 temp_Pos[1] = tPos[1];
1487
1488 /* Shift by pos -/+ on z */
1489 temp_Pos[2] -= dpos;
1490 azb1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1491 temp_Pos[2] = tPos[2];
1492
1493 temp_Pos[2] += dpos;
1494 azt1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1495 temp_Pos[2] = tPos[2];
1496
1497 /* Calculate the final value */
1498 dSA[0] = (axt1-axb1)/(2.0 * dpos);
1499 dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1500 dSA[2] = (azt1-azb1)/(2.0 * dpos);
1501}
1502
1503VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius) {
1504
1505 int i;
1506 int npts[3];
1507
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;
1512
1513 sav = 0.0;
1514 vol = 1.0;
1515 vol_density = 2.0;
1516
1517 lower_corner = clist->lower_corner;
1518 upper_corner = clist->upper_corner;
1519
1520 for (i=0; i<3; i++) {
1521 len = upper_corner[i] - lower_corner[i];
1522 vol *= len;
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) {
1527 if (apolparm->setgrid) {
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]);
1531 }
1532 spacs[i] = apolparm->grid[i];
1533
1534 }
1535 }
1536 }
1537
1538 for (x=lower_corner[0]; x<=upper_corner[0]; x=x+spacs[0]) {
1539 if ( VABS(x - lower_corner[0]) < VSMALL) {
1540 wx = 0.5;
1541 } else if ( VABS(x - upper_corner[0]) < VSMALL) {
1542 wx = 0.5;
1543 } else {
1544 wx = 1.0;
1545 }
1546 vec[0] = x;
1547 for (y=lower_corner[1]; y<=upper_corner[1]; y=y+spacs[1]) {
1548 if ( VABS(y - lower_corner[1]) < VSMALL) {
1549 wy = 0.5;
1550 } else if ( VABS(y - upper_corner[1]) < VSMALL) {
1551 wy = 0.5;
1552 } else {
1553 wy = 1.0;
1554 }
1555 vec[1] = y;
1556 for (z=lower_corner[2]; z<=upper_corner[2]; z=z+spacs[2]) {
1557 if ( VABS(z - lower_corner[2]) < VSMALL) {
1558 wz = 0.5;
1559 } else if ( VABS(z - upper_corner[2]) < VSMALL) {
1560 wz = 0.5;
1561 } else {
1562 wz = 1.0;
1563 }
1564 vec[2] = z;
1565
1566 w = wx*wy*wz;
1567
1568 sav += (w*(1.0-Vacc_ivdwAcc(thee, vec, radius)));
1569
1570 } /* z loop */
1571 } /* y loop */
1572 } /* x loop */
1573
1574 w = spacs[0]*spacs[1]*spacs[2];
1575 sav *= w;
1576
1577 return sav;
1578}
1579
1580int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist,
1581 Vclist *clist, int iatom, double *value) {
1582
1583 int i;
1584 int npts[3];
1585 int pad = 14;
1586
1587 int xmin, ymin, zmin;
1588 int xmax, ymax, zmax;
1589
1590 double sigma6, sigma12;
1591
1592 double spacs[3], vec[3];
1593 double w, wx, wy, wz, len, fn, x, y, z, vol;
1594 double x2,y2,z2,r;
1595 double vol_density, energy, rho, srad;
1596 double psig, epsilon, watepsilon, sigma, watsigma, eni, chi;
1597
1598 double *pos;
1599 double *lower_corner, *upper_corner;
1600
1601 Vatom *atom = VNULL;
1602 VASSERT(apolparm != VNULL);
1603
1604 energy = 0.0;
1605 vol = 1.0;
1606 vol_density = 2.0;
1607
1608 lower_corner = clist->lower_corner;
1609 upper_corner = clist->upper_corner;
1610
1611 atom = Valist_getAtom(alist, iatom);
1612 pos = Vatom_getPosition(atom);
1613
1614 /* Note: these are the original temporary water parameters... they have been
1615 replaced by entries in a parameter file:
1616 watsigma = 1.7683;
1617 watepsilon = 0.1521;
1618 watepsilon = watepsilon*4.184;
1619 */
1620
1621 srad = apolparm->srad;
1622 rho = apolparm->bconc;
1623 watsigma = apolparm->watsigma;
1624 watepsilon = apolparm->watepsilon;
1625 psig = atom->radius;
1626 epsilon = atom->epsilon;
1627 sigma = psig + watsigma;
1628 epsilon = VSQRT((epsilon * watepsilon));
1629
1630 /* parameters */
1631 sigma6 = VPOW(sigma,6);
1632 sigma12 = VPOW(sigma,12);
1633 /* OPLS-style radius: double sigmar = sigma*VPOW(2, (1.0/6.0)); */
1634
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;
1641
1642 for (i=0; i<3; i++) {
1643 len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1644 vol *= len;
1645 fn = len*vol_density + 1;
1646 npts[i] = (int)ceil(fn);
1647 spacs[i] = 0.5;
1648 if (apolparm->setgrid) {
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]);
1652 }
1653 spacs[i] = apolparm->grid[i];
1654 }
1655 }
1656
1657 for (x=xmin; x<=xmax; x=x+spacs[0]) {
1658 if ( VABS(x - xmin) < VSMALL) {
1659 wx = 0.5;
1660 } else if ( VABS(x - xmax) < VSMALL) {
1661 wx = 0.5;
1662 } else {
1663 wx = 1.0;
1664 }
1665 vec[0] = x;
1666 for (y=ymin; y<=ymax; y=y+spacs[1]) {
1667 if ( VABS(y - ymin) < VSMALL) {
1668 wy = 0.5;
1669 } else if ( VABS(y - ymax) < VSMALL) {
1670 wy = 0.5;
1671 } else {
1672 wy = 1.0;
1673 }
1674 vec[1] = y;
1675 for (z=zmin; z<=zmax; z=z+spacs[2]) {
1676 if ( VABS(z - zmin) < VSMALL) {
1677 wz = 0.5;
1678 } else if ( VABS(z - zmax) < VSMALL) {
1679 wz = 0.5;
1680 } else {
1681 wz = 1.0;
1682 }
1683 vec[2] = z;
1684
1685 w = wx*wy*wz;
1686
1687 chi = Vacc_ivdwAcc(thee, vec, srad);
1688
1689 if (VABS(chi) > VSMALL) {
1690
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);
1695
1696 if (r <= 14 && r >= sigma) {
1697 eni = chi*rho*epsilon*(-2.0*sigma6/VPOW(r,6)+sigma12/VPOW(r,12));
1698 }else if (r <= 14){
1699 eni = -1.0*epsilon*chi*rho;
1700 }else{
1701 eni = 0.0;
1702 }
1703 }else{
1704 eni = 0.0;
1705 }
1706
1707 energy += eni*w;
1708
1709 } /* z loop */
1710 } /* y loop */
1711 } /* x loop */
1712
1713 w = spacs[0]*spacs[1]*spacs[2];
1714 energy *= w;
1715
1716 *value = energy;
1717
1718 return VRC_SUCCESS;
1719}
1720
1721VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist,
1722 Vclist *clist){
1723
1724 int iatom;
1725 int rc = 0;
1726
1727 double energy = 0.0;
1728 double tenergy = 0.0;
1729 double rho = apolparm->bconc;
1730
1731 /* Do a sanity check to make sure that watepsilon and watsigma are set
1732 * If not, return with an error. */
1733 if(apolparm->setwat == 0){
1734 Vnm_print(2,"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1735 return VRC_FAILURE;
1736 }
1737
1738 if (VABS(rho) < VSMALL) {
1739 apolparm->wcaEnergy = tenergy;
1740 return 1;
1741 }
1742
1743 for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++){
1744 rc = Vacc_wcaEnergyAtom(acc, apolparm, alist, clist, iatom, &energy);
1745 if(rc == 0) return 0;
1746
1747 tenergy += energy;
1748 }
1749
1750 apolparm->wcaEnergy = tenergy;
1751
1752 return VRC_SUCCESS;
1753
1754}
1755
1756VPUBLIC int Vacc_wcaForceAtom(Vacc *thee,
1757 APOLparm *apolparm,
1758 Vclist *clist,
1759 Vatom *atom,
1760 double *force
1761 ){
1762 int i,
1763 si,
1764 npts[3],
1765 pad = 14,
1766 xmin,
1767 ymin,
1768 zmin,
1769 xmax,
1770 ymax,
1771 zmax;
1772
1773 double sigma6,
1774 sigma12,
1775 spacs[3],
1776 vec[3],
1777 fpt[3],
1778 w,
1779 wx,
1780 wy,
1781 wz,
1782 len,
1783 fn,
1784 x,
1785 y,
1786 z,
1787 vol,
1788 x2,
1789 y2,
1790 z2,
1791 r,
1792 vol_density,
1793 fo,
1794 rho,
1795 srad,
1796 psig,
1797 epsilon,
1798 watepsilon,
1799 sigma,
1800 watsigma,
1801 chi,
1802 *pos,
1803 *lower_corner,
1804 *upper_corner;
1805
1806 /* Allocate needed variables now that we've asserted required conditions. */
1807 time_t ts;
1808 ts = clock();
1809
1810 VASSERT(apolparm != VNULL);
1811
1812 /* Do a sanity check to make sure that watepsilon and watsigma are set
1813 * If not, return with an error. */
1814 if(apolparm->setwat == 0){
1815 Vnm_print(2,"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1816 return VRC_FAILURE;
1817 }
1818
1819 vol = 1.0;
1820 vol_density = 2.0;
1821
1822 lower_corner = clist->lower_corner;
1823 upper_corner = clist->upper_corner;
1824
1825 pos = Vatom_getPosition(atom);
1826
1827 srad = apolparm->srad;
1828 rho = apolparm->bconc;
1829 watsigma = apolparm->watsigma;
1830 watepsilon = apolparm->watepsilon;
1831
1832 psig = atom->radius;
1833 epsilon = atom->epsilon;
1834 sigma = psig + watsigma;
1835 epsilon = VSQRT((epsilon * watepsilon));
1836
1837 /* parameters */
1838 sigma6 = VPOW(sigma,6);
1839 sigma12 = VPOW(sigma,12);
1840 /* OPLS-style radius: double sigmar = sigma*VPOW(2, (1.0/6.0)); */
1841
1842 for (i=0; i<3; i++) {
1843 len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1844 vol *= len;
1845 fn = len*vol_density + 1;
1846 npts[i] = (int)ceil(fn);
1847 spacs[i] = 0.5;
1848 force[i] = 0.0;
1849 if (apolparm->setgrid) {
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]);
1853 }
1854 spacs[i] = apolparm->grid[i];
1855 }
1856 }
1857
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;
1864
1865 for (x=xmin; x<=xmax; x=x+spacs[0]) {
1866 if ( VABS(x - xmin) < VSMALL) {
1867 wx = 0.5;
1868 } else if ( VABS(x - xmax) < VSMALL) {
1869 wx = 0.5;
1870 } else {
1871 wx = 1.0;
1872 }
1873 vec[0] = x;
1874 for (y=ymin; y<=ymax; y=y+spacs[1]) {
1875 if ( VABS(y - ymin) < VSMALL) {
1876 wy = 0.5;
1877 } else if ( VABS(y - ymax) < VSMALL) {
1878 wy = 0.5;
1879 } else {
1880 wy = 1.0;
1881 }
1882 vec[1] = y;
1883 for (z=zmin; z<=zmax; z=z+spacs[2]) {
1884 if ( VABS(z - zmin) < VSMALL) {
1885 wz = 0.5;
1886 } else if ( VABS(z - zmax) < VSMALL) {
1887 wz = 0.5;
1888 } else {
1889 wz = 1.0;
1890 }
1891 vec[2] = z;
1892
1893 w = wx*wy*wz;
1894
1895 chi = Vacc_ivdwAcc(thee, vec, srad);
1896
1897 if (chi != 0.0) {
1898
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);
1903
1904 if (r <= 14 && r >= sigma){
1905
1906 fo = 12.0*chi*rho*epsilon*(sigma6/VPOW(r,7)-sigma12/VPOW(r,13));
1907
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;
1911
1912 }else {
1913 for (si=0; si < 3; si++) fpt[si] = 0.0;
1914 }
1915 }else {
1916 for (si=0; si < 3; si++) fpt[si] = 0.0;
1917 }
1918
1919 for(i=0;i<3;i++){
1920 force[i] += (w*fpt[i]);
1921 }
1922
1923 } /* z loop */
1924 } /* y loop */
1925 } /* x loop */
1926
1927 w = spacs[0]*spacs[1]*spacs[2];
1928 for(i=0;i<3;i++) force[i] *= w;
1929
1930 return VRC_SUCCESS;
1931}
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
Definition vacc.c:1721
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...
Definition vacc.c:316
VPUBLIC int Vacc_ctor2(Vacc *thee, Valist *alist, Vclist *clist, double surf_density)
FORTRAN stub to construct the accessibility object.
Definition vacc.c:213
int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist, Vclist *clist, int iatom, double *value)
Calculate the WCA energy for an atom.
Definition vacc.c:1580
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...
Definition vacc.c:377
VPUBLIC unsigned long int Vacc_memChk(Vacc *thee)
Get number of bytes in this object and its members.
Definition vacc.c:63
VPUBLIC void VaccSurf_dtor2(VaccSurf *thee)
Destroy the surface object.
Definition vacc.c:853
VPUBLIC void Vacc_totalAtomdSAV(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA, Vclist *clist)
Total solvent accessible volume.
Definition vacc.c:1448
VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility.
Definition vacc.c:608
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...
Definition vacc.c:1099
VPUBLIC void VaccSurf_dtor(VaccSurf **thee)
Destroy the surface object and free its memory.
Definition vacc.c:839
VPUBLIC void Vacc_totalAtomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Testing purposes only.
Definition vacc.c:1389
VPUBLIC void Vacc_splineAccGrad(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, double *grad)
Report gradient of spline-based accessibility.
Definition vacc.c:561
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
Definition vacc.c:245
VPUBLIC VaccSurf * VaccSurf_ctor(Vmem *mem, double probe_radius, int nsphere)
Allocate and construct the surface object; do not assign surface points to positions.
Definition vacc.c:803
VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad)
Report spline-based accessibility.
Definition vacc.c:528
VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
Definition vacc.c:774
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
Definition vacc.c:780
VPUBLIC VaccSurf * VaccSurf_refSphere(Vmem *mem, int npts)
Set up an array of points for a reference sphere of unit radius.
Definition vacc.c:926
VPUBLIC double Vacc_fastMolAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility quickly.
Definition vacc.c:637
VPUBLIC double Vacc_SASA(Vacc *thee, double radius)
Build the solvent accessible surface (SAS) and calculate the solvent accessible surface area.
Definition vacc.c:713
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...
Definition vacc.c:813
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...
Definition vacc.c:1006
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.
Definition vacc.c:868
VPUBLIC int Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
Definition vacc.c:1756
VPUBLIC VaccSurf * Vacc_atomSASPoints(Vacc *thee, double radius, Vatom *atom)
Get the set of points for this atom's solvent-accessible surface.
Definition vacc.c:982
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
Definition vacc.c:1503
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
Definition vacc.c:132
VPUBLIC void Vacc_dtor2(Vacc *thee)
FORTRAN stub to destroy object.
Definition vacc.c:255
VPUBLIC double Vacc_splineAccAtom(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, Vatom *atom)
Report spline-based accessibility for a given atom.
Definition vacc.c:438
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition valist.c:115
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition valist.c:105
VPUBLIC double Vatom_getAtomID(Vatom *thee)
Get atom ID.
Definition vatom.c:84
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
Definition vatom.c:105
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition vatom.c:63
VPUBLIC VclistCell * Vclist_getCell(Vclist *thee, double pos[VAPBS_DIM])
Return cell corresponding to specified position or return VNULL.
Definition vclist.c:423
VPUBLIC double Vclist_maxRadius(Vclist *thee)
Get the max probe radius value (in A) the cell list was constructed with.
Definition vclist.c:68
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition vhal.h:556
#define VAPBS_DIM
Our dimension.
Definition vhal.h:402
@ VRC_FAILURE
Definition vhal.h:69
@ VRC_SUCCESS
Definition vhal.h:70
#define Vunit_pi
Pi.
Definition vunit.h:104
Parameter structure for APOL-specific variables from input files.
Definition apolparm.h:129
double grid[3]
Definition apolparm.h:133
int setwat
Definition apolparm.h:180
double watepsilon
Definition apolparm.h:174
int setgrid
Definition apolparm.h:134
double bconc
Definition apolparm.h:139
double wcaEnergy
Definition apolparm.h:177
double watsigma
Definition apolparm.h:173
double srad
Definition apolparm.h:154
Oracle for solvent- and ion-accessibility around a biomolecule.
Definition vacc.h:108
Vmem * mem
Definition vacc.h:110
int * atomFlags
Definition vacc.h:113
double surf_density
Definition vacc.h:122
Vclist * clist
Definition vacc.h:112
Valist * alist
Definition vacc.h:111
VaccSurf * refSphere
Definition vacc.h:116
VaccSurf ** surf
Definition vacc.h:117
Surface object list of per-atom surface points.
Definition vacc.h:84
int npts
Definition vacc.h:92
Vmem * mem
Definition vacc.h:85
double * zpts
Definition vacc.h:88
double * ypts
Definition vacc.h:87
double probe_radius
Definition vacc.h:93
double * xpts
Definition vacc.h:86
double area
Definition vacc.h:91
char * bpts
Definition vacc.h:89
Container class for list of atom objects.
Definition valist.h:78
Contains public data members for Vatom class/module.
Definition vatom.h:84
double radius
Definition vatom.h:87
double epsilon
Definition vatom.h:91
int id
Definition vatom.h:93
double position[3]
Definition vatom.h:86
Atom cell list cell.
Definition vclist.h:101
Vatom ** atoms
Definition vclist.h:102
int natoms
Definition vclist.h:103
Atom cell list.
Definition vclist.h:117
double upper_corner[VAPBS_DIM]
Definition vclist.h:127
double lower_corner[VAPBS_DIM]
Definition vclist.h:126
VPRIVATE double splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad, VclistCell *cell)
Fast spline-based surface computation subroutine.
Definition vacc.c:493
VPRIVATE int Vacc_storeParms(Vacc *thee, Valist *alist, Vclist *clist, double surf_density)
Definition vacc.c:148
VPRIVATE int Vacc_allocate(Vacc *thee)
Definition vacc.c:193
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...
Definition vacc.c:80
Contains declarations for class Vacc.