APBS  1.4.1
vacc.c
Go to the documentation of this file.
1 
57 #include "vacc.h"
58 
59 VEMBED(rcsid="$Id$")
60 
61 #if !defined(VINLINE_VACC)
62 
63 VPUBLIC 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 
80 VPRIVATE 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 
132 VPUBLIC 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 
148 VPRIVATE 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 
193 VPRIVATE 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 
213 VPUBLIC 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 
245 VPUBLIC 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 
255 VPUBLIC 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 
277 VPUBLIC 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 
307 VPUBLIC 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 
316 VPUBLIC void Vacc_splineAccGradAtomNorm(Vacc *thee,
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 
438 VPUBLIC 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 
493 VPRIVATE 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 
528 VPUBLIC 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 
561 VPUBLIC 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 
608 VPUBLIC 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 
637 VPUBLIC 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)
679 VPUBLIC 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 
713 VPUBLIC 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 
774 VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius) {
775 
776  return Vacc_SASA(thee, radius);
777 
778 }
779 
780 VPUBLIC 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 
803 VPUBLIC VaccSurf* VaccSurf_ctor(Vmem *mem, double probe_radius, int nsphere) {
804  VaccSurf *thee;
805 
806  //thee = Vmem_malloc(mem, 1, sizeof(Vacc) );
807  if (nsphere >= MAX_SPHERE_PTS) {
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");
810  VASSERT(0);
811  }
812  thee = (VaccSurf*)calloc(1,sizeof(Vacc));
813  VASSERT( VaccSurf_ctor2(thee, mem, probe_radius, nsphere) );
814 
815  return thee;
816 }
817 
818 VPUBLIC int VaccSurf_ctor2(VaccSurf *thee, Vmem *mem, double probe_radius,
819  int nsphere) {
820 
821  if (thee == VNULL)
822  return 0;
823 
824  thee->mem = mem;
825  thee->npts = nsphere;
826  thee->probe_radius = probe_radius;
827  thee->area = 0.0;
828 
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));
834  } else {
835  thee->xpts = VNULL;
836  thee->ypts = VNULL;
837  thee->zpts = VNULL;
838  thee->bpts = VNULL;
839  }
840 
841  return 1;
842 }
843 
844 VPUBLIC void VaccSurf_dtor(VaccSurf **thee) {
845 
846  Vmem *mem;
847 
848  if ((*thee) != VNULL) {
849  mem = (*thee)->mem;
850  VaccSurf_dtor2(*thee);
851  //Vmem_free(mem, 1, sizeof(VaccSurf), (void **)thee);
852  free(*thee);
853  (*thee) = VNULL;
854  }
855 
856 }
857 
858 VPUBLIC void VaccSurf_dtor2(VaccSurf *thee) {
859 
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));
869  }
870 
871 }
872 
873 VPUBLIC VaccSurf* Vacc_atomSurf(Vacc *thee,
874  Vatom *atom,
875  VaccSurf *ref,
876  double prad
877  ) {
878 
879  VaccSurf *surf;
880  int i,
881  j,
882  npts,
883  atomID;
884  double arad,
885  rad,
886  pos[3],
887  *apos;
888  char bpts[MAX_SPHERE_PTS];
889 
890  /* Get atom information */
891  arad = Vatom_getRadius(atom);
892  apos = Vatom_getPosition(atom);
893  atomID = Vatom_getAtomID(atom);
894 
895  if (arad < VSMALL) {
896  return VaccSurf_ctor(thee->mem, prad, 0);
897  }
898 
899  rad = arad + prad;
900 
901  /* Determine which points will contribute */
902  npts = 0;
903  for (i=0; i<ref->npts; i++) {
904  /* Reset point flag: zero-radius atoms do not contribute */
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];
908  if (ivdwAccExclus(thee, pos, prad, atomID)) {
909  npts++;
910  bpts[i] = 1;
911  } else {
912  bpts[i] = 0;
913  }
914  }
915 
916  /* Allocate space for the points */
917  surf = VaccSurf_ctor(thee->mem, prad, npts);
918 
919  /* Assign the points */
920  j = 0;
921  for (i=0; i<ref->npts; i++) {
922  if (bpts[i]) {
923  surf->bpts[j] = 1;
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];
927  j++;
928  }
929  }
930 
931  /* Assign the area */
932  surf->area = 4.0*VPI*rad*rad*((double)(surf->npts))/((double)(ref->npts));
933 
934  return surf;
935 
936 }
937 
938 VPUBLIC VaccSurf* VaccSurf_refSphere(Vmem *mem, int npts) {
939 
940  VaccSurf *surf;
941  int nactual, i, itheta, ntheta, iphi, nphimax, nphi;
942  double frac;
943  double sintheta, costheta, theta, dtheta;
944  double sinphi, cosphi, phi, dphi;
945 
946  /* Setup "constants" */
947  frac = ((double)(npts))/4.0;
948  ntheta = VRINT(VSQRT(Vunit_pi*frac));
949  dtheta = Vunit_pi/((double)(ntheta));
950  nphimax = 2*ntheta;
951 
952  /* Count the actual number of points to be used */
953  nactual = 0;
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);
959  nactual += nphi;
960  }
961 
962  /* Allocate space for the points */
963  surf = VaccSurf_ctor(mem, 1.0, nactual);
964 
965  /* Clear out the boolean array */
966  for (i=0; i<nactual; i++) surf->bpts[i] = 1;
967 
968  /* Assign the points */
969  nactual = 0;
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);
975  if (nphi != 0) {
976  dphi = 2*Vunit_pi/((double)(nphi));
977  for (iphi=0; iphi<nphi; iphi++) {
978  phi = dphi*((double)(iphi));
979  sinphi = VSIN(phi);
980  cosphi = VCOS(phi);
981  surf->xpts[nactual] = cosphi * sintheta;
982  surf->ypts[nactual] = sinphi * sintheta;
983  surf->zpts[nactual] = costheta;
984  nactual++;
985  }
986  }
987  }
988 
989  surf->npts = nactual;
990 
991  return surf;
992 }
993 
994 VPUBLIC VaccSurf* Vacc_atomSASPoints(Vacc *thee, double radius,
995  Vatom *atom) {
996 
997  VaccSurf *asurf = VNULL;
998  int id;
999 
1000  if (thee->surf == VNULL) Vacc_SASA(thee, radius);
1001  id = Vatom_getAtomID(atom);
1002 
1003  asurf = thee->surf[id];
1004 
1005  /* See if this surface needs to be rebuilt */
1006  if (asurf->probe_radius != radius) {
1007  Vnm_print(2, "Vacc_SASA: Warning -- probe radius changed from %g to %g!\n",
1008  asurf->probe_radius, radius);
1009  VaccSurf_dtor2(asurf);
1010  thee->surf[id] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
1011  asurf = thee->surf[id];
1012  }
1013 
1014  return asurf;
1015 
1016 }
1017 
1018 VPUBLIC void Vacc_splineAccGradAtomNorm4(Vacc *thee, double center[VAPBS_DIM],
1019  double win, double infrad, Vatom *atom, double *grad) {
1020 
1021  int i;
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;
1027  double mychi = 1.0; /* Char. func. value for given atom */
1028 
1029  VASSERT(thee != NULL);
1030 
1031  /* The grad is zero by default */
1032  for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
1033 
1034  /* *** CALCULATE THE CHARACTERISTIC FUNCTION VALUE FOR THIS ATOM AND THE
1035  * *** MAGNITUDE OF THE FORCE *** */
1036  apos = Vatom_getPosition(atom);
1037  /* Zero-radius atoms don't contribute */
1038  if (Vatom_getRadius(atom) > 0.0) {
1039 
1040  arad = Vatom_getRadius(atom);
1041  arad = arad + infrad;
1042  b = arad - win;
1043  e = arad + win;
1044 
1045  e2 = e * e;
1046  e3 = e2 * e;
1047  e4 = e3 * e;
1048  e5 = e4 * e;
1049  e6 = e5 * e;
1050  e7 = e6 * e;
1051  b2 = b * b;
1052  b3 = b2 * b;
1053  b4 = b3 * b;
1054  b5 = b4 * b;
1055  b6 = b5 * b;
1056  b7 = b6 * b;
1057 
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;
1067  c7 = -20.0/denom;
1068 
1069  dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1070  + VSQR(apos[2]-center[2]));
1071 
1072  /* If we're inside an atom, the entire characteristic function
1073  * will be zero and the grad will be zero, so we can stop */
1074  if (dist < (arad - win)) return;
1075  /* Likewise, if we're outside the smoothing window, the characteristic
1076  * function is unity and the grad will be zero, so we can stop */
1077  else if (dist > (arad + win)) return;
1078  /* Account for floating point error at the border
1079  * NAB: COULDN'T THESE TESTS BE COMBINED AS BELOW
1080  * (Vacc_splineAccAtom)? */
1081  else if ((VABS(dist - (arad - win)) < VSMALL) ||
1082  (VABS(dist - (arad + win)) < VSMALL)) return;
1083  /* If we're inside the smoothing window */
1084  else {
1085  sm = dist;
1086  sm2 = sm * sm;
1087  sm3 = sm2 * sm;
1088  sm4 = sm3 * sm;
1089  sm5 = sm4 * sm;
1090  sm6 = sm5 * sm;
1091  sm7 = sm6 * sm;
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;
1096  if (mychi <= 0.0) {
1097  /* Avoid numerical round off errors */
1098  return;
1099  } else if (mychi > 1.0) {
1100  /* Avoid numerical round off errors */
1101  mychi = 1.0;
1102  }
1103  }
1104  /* Now assemble the grad vector */
1105  VASSERT(mychi > 0.0);
1106  for (i=0; i<VAPBS_DIM; i++)
1107  grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1108  }
1109 }
1110 
1111 VPUBLIC void Vacc_splineAccGradAtomNorm3(Vacc *thee, double center[VAPBS_DIM],
1112  double win, double infrad, Vatom *atom, double *grad) {
1113 
1114  int i;
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;
1120  double mychi = 1.0; /* Char. func. value for given atom */
1121 
1122  VASSERT(thee != NULL);
1123 
1124  /* The grad is zero by default */
1125  for (i=0; i<VAPBS_DIM; i++) grad[i] = 0.0;
1126 
1127  /* *** CALCULATE THE CHARACTERISTIC FUNCTION VALUE FOR THIS ATOM AND THE
1128  * *** MAGNITUDE OF THE FORCE *** */
1129  apos = Vatom_getPosition(atom);
1130  /* Zero-radius atoms don't contribute */
1131  if (Vatom_getRadius(atom) > 0.0) {
1132 
1133  arad = Vatom_getRadius(atom);
1134  arad = arad + infrad;
1135  b = arad - win;
1136  e = arad + win;
1137 
1138  e2 = e * e;
1139  e3 = e2 * e;
1140  e4 = e3 * e;
1141  e5 = e4 * e;
1142  b2 = b * b;
1143  b3 = b2 * b;
1144  b4 = b3 * b;
1145  b5 = b4 * b;
1146 
1147  denom = pow((e - b), 5.0);
1148  c0 = -10.0*e2*b3 + 5.0*e*b4 - b5;
1149  c1 = 30.0*e2*b2;
1150  c2 = -30.0*(e2*b + e*b2);
1151  c3 = 10.0*(e2 + 4.0*e*b + b2);
1152  c4 = -15.0*(e + b);
1153  c5 = 6;
1154  c0 = c0/denom;
1155  c1 = c1/denom;
1156  c2 = c2/denom;
1157  c3 = c3/denom;
1158  c4 = c4/denom;
1159  c5 = c5/denom;
1160 
1161  dist = VSQRT(VSQR(apos[0]-center[0]) + VSQR(apos[1]-center[1])
1162  + VSQR(apos[2]-center[2]));
1163 
1164  /* If we're inside an atom, the entire characteristic function
1165  * will be zero and the grad will be zero, so we can stop */
1166  if (dist < (arad - win)) return;
1167  /* Likewise, if we're outside the smoothing window, the characteristic
1168  * function is unity and the grad will be zero, so we can stop */
1169  else if (dist > (arad + win)) return;
1170  /* Account for floating point error at the border
1171  * NAB: COULDN'T THESE TESTS BE COMBINED AS BELOW
1172  * (Vacc_splineAccAtom)? */
1173  else if ((VABS(dist - (arad - win)) < VSMALL) ||
1174  (VABS(dist - (arad + win)) < VSMALL)) return;
1175  /* If we're inside the smoothing window */
1176  else {
1177  sm = dist;
1178  sm2 = sm * sm;
1179  sm3 = sm2 * sm;
1180  sm4 = sm3 * sm;
1181  sm5 = sm4 * sm;
1182  mychi = c0 + c1*sm + c2*sm2 + c3*sm3
1183  + c4*sm4 + c5*sm5;
1184  mygrad = c1 + 2.0*c2*sm + 3.0*c3*sm2 + 4.0*c4*sm3
1185  + 5.0*c5*sm4;
1186  if (mychi <= 0.0) {
1187  /* Avoid numerical round off errors */
1188  return;
1189  } else if (mychi > 1.0) {
1190  /* Avoid numerical round off errors */
1191  mychi = 1.0;
1192  }
1193  }
1194  /* Now assemble the grad vector */
1195  VASSERT(mychi > 0.0);
1196  for (i=0; i<VAPBS_DIM; i++)
1197  grad[i] = -(mygrad/mychi)*((center[i] - apos[i])/dist);
1198  }
1199 }
1200 
1201 /* ///////////////////////////////////////////////////////////////////////////
1202  // Routine: Vacc_atomdSAV
1203  //
1204  // Purpose: Calculates the vector valued atomic derivative of volume
1205  //
1206  // Args: radius The radius of the solvent probe in Angstroms
1207  // iatom Index of the atom in thee->alist
1208  //
1209  // Author: Jason Wagoner
1210  // Nathan Baker (original FORTRAN routine from UHBD by Brock Luty)
1212 VPUBLIC void Vacc_atomdSAV(Vacc *thee,
1213  double srad,
1214  Vatom *atom,
1215  double *dSA
1216  ) {
1217 
1218  int ipt, iatom;
1219 
1220  double area;
1221  double *tPos, tRad, vec[3];
1222  double dx,dy,dz;
1223  VaccSurf *ref;
1224  dx = 0.0;
1225  dy = 0.0;
1226  dz = 0.0;
1227  /* Get the atom information */
1228  ref = thee->refSphere;
1229  iatom = Vatom_getAtomID(atom);
1230 
1231  dSA[0] = 0.0;
1232  dSA[1] = 0.0;
1233  dSA[2] = 0.0;
1234 
1235  tPos = Vatom_getPosition(atom);
1236  tRad = Vatom_getRadius(atom);
1237 
1238  if(tRad == 0.0) return;
1239 
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];
1245  if (ivdwAccExclus(thee, vec, srad, iatom)) {
1246  dx = dx+vec[0]-tPos[0];
1247  dy = dy+vec[1]-tPos[1];
1248  dz = dz+vec[2]-tPos[2];
1249  }
1250  }
1251 
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);
1256  }
1257 
1258 }
1259 
1260 /* Note: This is purely test code to make certain that the dSASA code is
1261  behaving properly. This function should NEVER be called by anyone
1262  other than an APBS developer at Wash U.
1263 */
1264 VPRIVATE double Vacc_SASAPos(Vacc *thee, double radius) {
1265 
1266  int i, natom;
1267  double area;
1268  Vatom *atom;
1269  VaccSurf *asurf;
1270 
1271  natom = Valist_getNumberAtoms(thee->alist);
1272 
1273  /* Calculate the area */
1274  area = 0.0;
1275  for (i=0; i<natom; i++) {
1276  atom = Valist_getAtom(thee->alist, i);
1277  asurf = thee->surf[i];
1278 
1279  VaccSurf_dtor2(asurf);
1280  thee->surf[i] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
1281  asurf = thee->surf[i];
1282  area += (asurf->area);
1283  }
1284 
1285  return area;
1286 
1287 }
1288 
1289 VPRIVATE double Vacc_atomSASAPos(Vacc *thee,
1290  double radius,
1291  Vatom *atom, /* The atom being manipulated */
1292  int mode
1293  ) {
1294 
1295  VaccSurf *asurf;
1296  int id;
1297  static int warned = 0;
1298 
1299  if ((thee->surf == VNULL) || (mode == 1)){
1300  if(!warned){
1301  Vnm_print(2, "WARNING: Recalculating entire surface!!!!\n");
1302  warned = 1;
1303  }
1304  Vacc_SASAPos(thee, radius); // reinitialize before we can do anything about doing a calculation on a repositioned atom
1305  }
1306 
1307  id = Vatom_getAtomID(atom);
1308  asurf = thee->surf[id];
1309 
1310  VaccSurf_dtor(&asurf);
1311  thee->surf[id] = Vacc_atomSurf(thee, atom, thee->refSphere, radius);
1312  asurf = thee->surf[id];
1313 
1314  //printf("%s: Time elapsed: %f\n", __func__, ((double)clock() - ts) / CLOCKS_PER_SEC);
1315 
1316  return asurf->area;
1317 }
1318 
1319 /* ///////////////////////////////////////////////////////////////////////////
1320  // Routine: Vacc_atomdSASA
1321  //
1322  // Purpose: Calculates the derivative of surface area with respect to atomic
1323  // displacement using finite difference methods.
1324  //
1325  // Args: radius The radius of the solvent probe in Angstroms
1326  // iatom Index of the atom in thee->alist
1327  //
1328  // Author: Jason Wagoner
1329  // David Gohara
1330  // Nathan Baker (original FORTRAN routine from UHBD by Brock Luty)
1332 VPUBLIC void Vacc_atomdSASA(Vacc *thee,
1333  double dpos,
1334  double srad,
1335  Vatom *atom,
1336  double *dSA
1337  ) {
1338 
1339  double *temp_Pos,
1340  tPos[3],
1341  axb1,
1342  axt1,
1343  ayb1,
1344  ayt1,
1345  azb1,
1346  azt1;
1347  VaccSurf *ref;
1348 
1349  //printf("%s: entering\n", __func__);
1350  time_t ts;
1351  ts = clock();
1352 
1353  /* Get the atom information */
1354  ref = thee->refSphere;
1355  temp_Pos = Vatom_getPosition(atom); // Get a pointer to the position object. You actually manipulate the atom doing this...
1356 
1357  tPos[0] = temp_Pos[0];
1358  tPos[1] = temp_Pos[1];
1359  tPos[2] = temp_Pos[2];
1360 
1361  /* Shift by pos -/+ on x */
1362  temp_Pos[0] -= dpos;
1363  axb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1364  temp_Pos[0] = tPos[0];
1365 
1366  temp_Pos[0] += dpos;
1367  axt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1368  temp_Pos[0] = tPos[0];
1369 
1370  /* Shift by pos -/+ on y */
1371  temp_Pos[1] -= dpos;
1372  ayb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1373  temp_Pos[1] = tPos[1];
1374 
1375  temp_Pos[1] += dpos;
1376  ayt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1377  temp_Pos[1] = tPos[1];
1378 
1379  /* Shift by pos -/+ on z */
1380  temp_Pos[2] -= dpos;
1381  azb1 = Vacc_atomSASAPos(thee, srad, atom,0);
1382  temp_Pos[2] = tPos[2];
1383 
1384  temp_Pos[2] += dpos;
1385  azt1 = Vacc_atomSASAPos(thee, srad, atom,0);
1386  temp_Pos[2] = tPos[2];
1387 
1388  /* Reset the atom SASA to zero displacement */
1389  Vacc_atomSASAPos(thee, srad, atom,0);
1390 
1391  /* Calculate the final value */
1392  dSA[0] = (axt1-axb1)/(2.0 * dpos);
1393  dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1394  dSA[2] = (azt1-azb1)/(2.0 * dpos);
1395 }
1396 
1397 /* Note: This is purely test code to make certain that the dSASA code is
1398  behaving properly. This function should NEVER be called by anyone
1399  other than an APBS developer at Wash U.
1400 */
1401 VPUBLIC void Vacc_totalAtomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA) {
1402 
1403  int iatom;
1404  double *temp_Pos, tRad;
1405  double tPos[3];
1406  double axb1,axt1,ayb1,ayt1,azb1,azt1;
1407  VaccSurf *ref;
1408 
1409  /* Get the atom information */
1410  ref = thee->refSphere;
1411  temp_Pos = Vatom_getPosition(atom);
1412  tRad = Vatom_getRadius(atom);
1413  iatom = Vatom_getAtomID(atom);
1414 
1415  dSA[0] = 0.0;
1416  dSA[1] = 0.0;
1417  dSA[2] = 0.0;
1418 
1419  tPos[0] = temp_Pos[0];
1420  tPos[1] = temp_Pos[1];
1421  tPos[2] = temp_Pos[2];
1422 
1423  /* Shift by pos -/+ on x */
1424  temp_Pos[0] -= dpos;
1425  axb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1426  temp_Pos[0] = tPos[0];
1427 
1428  temp_Pos[0] += dpos;
1429  axt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1430  temp_Pos[0] = tPos[0];
1431 
1432  /* Shift by pos -/+ on y */
1433  temp_Pos[1] -= dpos;
1434  ayb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1435  temp_Pos[1] = tPos[1];
1436 
1437  temp_Pos[1] += dpos;
1438  ayt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1439  temp_Pos[1] = tPos[1];
1440 
1441  /* Shift by pos -/+ on z */
1442  temp_Pos[2] -= dpos;
1443  azb1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1444  temp_Pos[2] = tPos[2];
1445 
1446  temp_Pos[2] += dpos;
1447  azt1 = Vacc_atomSASAPos(thee, srad, atom, 1);
1448  temp_Pos[2] = tPos[2];
1449 
1450  /* Calculate the final value */
1451  dSA[0] = (axt1-axb1)/(2.0 * dpos);
1452  dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1453  dSA[2] = (azt1-azb1)/(2.0 * dpos);
1454 }
1455 
1456 /* Note: This is purely test code to make certain that the dSASA code is
1457  behaving properly. This function should NEVER be called by anyone
1458  other than an APBS developer at Wash U.
1459 */
1460 VPUBLIC void Vacc_totalAtomdSAV(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA, Vclist *clist) {
1461 
1462  int iatom;
1463  double *temp_Pos, tRad;
1464  double tPos[3];
1465  double axb1,axt1,ayb1,ayt1,azb1,azt1;
1466  VaccSurf *ref;
1467 
1468  /* Get the atom information */
1469  ref = thee->refSphere;
1470  temp_Pos = Vatom_getPosition(atom);
1471  tRad = Vatom_getRadius(atom);
1472  iatom = Vatom_getAtomID(atom);
1473 
1474  dSA[0] = 0.0;
1475  dSA[1] = 0.0;
1476  dSA[2] = 0.0;
1477 
1478  tPos[0] = temp_Pos[0];
1479  tPos[1] = temp_Pos[1];
1480  tPos[2] = temp_Pos[2];
1481 
1482  /* Shift by pos -/+ on x */
1483  temp_Pos[0] -= dpos;
1484  axb1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1485  temp_Pos[0] = tPos[0];
1486 
1487  temp_Pos[0] += dpos;
1488  axt1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1489  temp_Pos[0] = tPos[0];
1490 
1491  /* Shift by pos -/+ on y */
1492  temp_Pos[1] -= dpos;
1493  ayb1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1494  temp_Pos[1] = tPos[1];
1495 
1496  temp_Pos[1] += dpos;
1497  ayt1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1498  temp_Pos[1] = tPos[1];
1499 
1500  /* Shift by pos -/+ on z */
1501  temp_Pos[2] -= dpos;
1502  azb1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1503  temp_Pos[2] = tPos[2];
1504 
1505  temp_Pos[2] += dpos;
1506  azt1 = Vacc_totalSAV(thee,clist, VNULL, srad);
1507  temp_Pos[2] = tPos[2];
1508 
1509  /* Calculate the final value */
1510  dSA[0] = (axt1-axb1)/(2.0 * dpos);
1511  dSA[1] = (ayt1-ayb1)/(2.0 * dpos);
1512  dSA[2] = (azt1-azb1)/(2.0 * dpos);
1513 }
1514 
1515 VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius) {
1516 
1517  int i;
1518  int npts[3];
1519 
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;
1524 
1525  sav = 0.0;
1526  vol = 1.0;
1527  vol_density = 2.0;
1528 
1529  lower_corner = clist->lower_corner;
1530  upper_corner = clist->upper_corner;
1531 
1532  for (i=0; i<3; i++) {
1533  len = upper_corner[i] - lower_corner[i];
1534  vol *= len;
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) {
1539  if (apolparm->setgrid) {
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]);
1543  }
1544  spacs[i] = apolparm->grid[i];
1545 
1546  }
1547  }
1548  }
1549 
1550  for (x=lower_corner[0]; x<=upper_corner[0]; x=x+spacs[0]) {
1551  if ( VABS(x - lower_corner[0]) < VSMALL) {
1552  wx = 0.5;
1553  } else if ( VABS(x - upper_corner[0]) < VSMALL) {
1554  wx = 0.5;
1555  } else {
1556  wx = 1.0;
1557  }
1558  vec[0] = x;
1559  for (y=lower_corner[1]; y<=upper_corner[1]; y=y+spacs[1]) {
1560  if ( VABS(y - lower_corner[1]) < VSMALL) {
1561  wy = 0.5;
1562  } else if ( VABS(y - upper_corner[1]) < VSMALL) {
1563  wy = 0.5;
1564  } else {
1565  wy = 1.0;
1566  }
1567  vec[1] = y;
1568  for (z=lower_corner[2]; z<=upper_corner[2]; z=z+spacs[2]) {
1569  if ( VABS(z - lower_corner[2]) < VSMALL) {
1570  wz = 0.5;
1571  } else if ( VABS(z - upper_corner[2]) < VSMALL) {
1572  wz = 0.5;
1573  } else {
1574  wz = 1.0;
1575  }
1576  vec[2] = z;
1577 
1578  w = wx*wy*wz;
1579 
1580  sav += (w*(1.0-Vacc_ivdwAcc(thee, vec, radius)));
1581 
1582  } /* z loop */
1583  } /* y loop */
1584  } /* x loop */
1585 
1586  w = spacs[0]*spacs[1]*spacs[2];
1587  sav *= w;
1588 
1589  return sav;
1590 }
1591 
1592 int Vacc_wcaEnergyAtom(Vacc *thee, APOLparm *apolparm, Valist *alist,
1593  Vclist *clist, int iatom, double *value) {
1594 
1595  int i;
1596  int npts[3];
1597  int pad = 14;
1598 
1599  int xmin, ymin, zmin;
1600  int xmax, ymax, zmax;
1601 
1602  double sigma6, sigma12;
1603 
1604  double spacs[3], vec[3];
1605  double w, wx, wy, wz, len, fn, x, y, z, vol;
1606  double x2,y2,z2,r;
1607  double vol_density, energy, rho, srad;
1608  double psig, epsilon, watepsilon, sigma, watsigma, eni, chi;
1609 
1610  double *pos;
1611  double *lower_corner, *upper_corner;
1612 
1613  Vatom *atom = VNULL;
1614  VASSERT(apolparm != VNULL);
1615 
1616  energy = 0.0;
1617  vol = 1.0;
1618  vol_density = 2.0;
1619 
1620  lower_corner = clist->lower_corner;
1621  upper_corner = clist->upper_corner;
1622 
1623  atom = Valist_getAtom(alist, iatom);
1624  pos = Vatom_getPosition(atom);
1625 
1626  /* Note: these are the original temporary water parameters... they have been
1627  replaced by entries in a parameter file:
1628  watsigma = 1.7683;
1629  watepsilon = 0.1521;
1630  watepsilon = watepsilon*4.184;
1631  */
1632 
1633  srad = apolparm->srad;
1634  rho = apolparm->bconc;
1635  watsigma = apolparm->watsigma;
1636  watepsilon = apolparm->watepsilon;
1637  psig = atom->radius;
1638  epsilon = atom->epsilon;
1639  sigma = psig + watsigma;
1640  epsilon = VSQRT((epsilon * watepsilon));
1641 
1642  /* parameters */
1643  sigma6 = VPOW(sigma,6);
1644  sigma12 = VPOW(sigma,12);
1645  /* OPLS-style radius: double sigmar = sigma*VPOW(2, (1.0/6.0)); */
1646 
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;
1653 
1654  for (i=0; i<3; i++) {
1655  len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1656  vol *= len;
1657  fn = len*vol_density + 1;
1658  npts[i] = (int)ceil(fn);
1659  spacs[i] = 0.5;
1660  if (apolparm->setgrid) {
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]);
1664  }
1665  spacs[i] = apolparm->grid[i];
1666  }
1667  }
1668 
1669  for (x=xmin; x<=xmax; x=x+spacs[0]) {
1670  if ( VABS(x - xmin) < VSMALL) {
1671  wx = 0.5;
1672  } else if ( VABS(x - xmax) < VSMALL) {
1673  wx = 0.5;
1674  } else {
1675  wx = 1.0;
1676  }
1677  vec[0] = x;
1678  for (y=ymin; y<=ymax; y=y+spacs[1]) {
1679  if ( VABS(y - ymin) < VSMALL) {
1680  wy = 0.5;
1681  } else if ( VABS(y - ymax) < VSMALL) {
1682  wy = 0.5;
1683  } else {
1684  wy = 1.0;
1685  }
1686  vec[1] = y;
1687  for (z=zmin; z<=zmax; z=z+spacs[2]) {
1688  if ( VABS(z - zmin) < VSMALL) {
1689  wz = 0.5;
1690  } else if ( VABS(z - zmax) < VSMALL) {
1691  wz = 0.5;
1692  } else {
1693  wz = 1.0;
1694  }
1695  vec[2] = z;
1696 
1697  w = wx*wy*wz;
1698 
1699  chi = Vacc_ivdwAcc(thee, vec, srad);
1700 
1701  if (VABS(chi) > VSMALL) {
1702 
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);
1707 
1708  if (r <= 14 && r >= sigma) {
1709  eni = chi*rho*epsilon*(-2.0*sigma6/VPOW(r,6)+sigma12/VPOW(r,12));
1710  }else if (r <= 14){
1711  eni = -1.0*epsilon*chi*rho;
1712  }else{
1713  eni = 0.0;
1714  }
1715  }else{
1716  eni = 0.0;
1717  }
1718 
1719  energy += eni*w;
1720 
1721  } /* z loop */
1722  } /* y loop */
1723  } /* x loop */
1724 
1725  w = spacs[0]*spacs[1]*spacs[2];
1726  energy *= w;
1727 
1728  *value = energy;
1729 
1730  return VRC_SUCCESS;
1731 }
1732 
1733 VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist,
1734  Vclist *clist){
1735 
1736  int iatom;
1737  int rc = 0;
1738 
1739  double energy = 0.0;
1740  double tenergy = 0.0;
1741  double rho = apolparm->bconc;
1742 
1743  /* Do a sanity check to make sure that watepsilon and watsigma are set
1744  * If not, return with an error. */
1745  if(apolparm->setwat == 0){
1746  Vnm_print(2,"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1747  return VRC_FAILURE;
1748  }
1749 
1750  if (VABS(rho) < VSMALL) {
1751  apolparm->wcaEnergy = tenergy;
1752  return 1;
1753  }
1754 
1755  for (iatom=0; iatom<Valist_getNumberAtoms(alist); iatom++){
1756  rc = Vacc_wcaEnergyAtom(acc, apolparm, alist, clist, iatom, &energy);
1757  if(rc == 0) return 0;
1758 
1759  tenergy += energy;
1760  }
1761 
1762  apolparm->wcaEnergy = tenergy;
1763 
1764  return VRC_SUCCESS;
1765 
1766 }
1767 
1768 VPUBLIC int Vacc_wcaForceAtom(Vacc *thee,
1769  APOLparm *apolparm,
1770  Vclist *clist,
1771  Vatom *atom,
1772  double *force
1773  ){
1774  int i,
1775  si,
1776  npts[3],
1777  pad = 14,
1778  xmin,
1779  ymin,
1780  zmin,
1781  xmax,
1782  ymax,
1783  zmax;
1784 
1785  double sigma6,
1786  sigma12,
1787  spacs[3],
1788  vec[3],
1789  fpt[3],
1790  w,
1791  wx,
1792  wy,
1793  wz,
1794  len,
1795  fn,
1796  x,
1797  y,
1798  z,
1799  vol,
1800  x2,
1801  y2,
1802  z2,
1803  r,
1804  vol_density,
1805  fo,
1806  rho,
1807  srad,
1808  psig,
1809  epsilon,
1810  watepsilon,
1811  sigma,
1812  watsigma,
1813  chi,
1814  *pos,
1815  *lower_corner,
1816  *upper_corner;
1817 
1818  /* Allocate needed variables now that we've asserted required conditions. */
1819  time_t ts;
1820  ts = clock();
1821 
1822  VASSERT(apolparm != VNULL);
1823 
1824  /* Do a sanity check to make sure that watepsilon and watsigma are set
1825  * If not, return with an error. */
1826  if(apolparm->setwat == 0){
1827  Vnm_print(2,"Vacc_wcaEnergy: Error. No value was set for watsigma and watepsilon.\n");
1828  return VRC_FAILURE;
1829  }
1830 
1831  vol = 1.0;
1832  vol_density = 2.0;
1833 
1834  lower_corner = clist->lower_corner;
1835  upper_corner = clist->upper_corner;
1836 
1837  pos = Vatom_getPosition(atom);
1838 
1839  srad = apolparm->srad;
1840  rho = apolparm->bconc;
1841  watsigma = apolparm->watsigma;
1842  watepsilon = apolparm->watepsilon;
1843 
1844  psig = atom->radius;
1845  epsilon = atom->epsilon;
1846  sigma = psig + watsigma;
1847  epsilon = VSQRT((epsilon * watepsilon));
1848 
1849  /* parameters */
1850  sigma6 = VPOW(sigma,6);
1851  sigma12 = VPOW(sigma,12);
1852  /* OPLS-style radius: double sigmar = sigma*VPOW(2, (1.0/6.0)); */
1853 
1854  for (i=0; i<3; i++) {
1855  len = (upper_corner[i] + pad) - (lower_corner[i] - pad);
1856  vol *= len;
1857  fn = len*vol_density + 1;
1858  npts[i] = (int)ceil(fn);
1859  spacs[i] = 0.5;
1860  force[i] = 0.0;
1861  if (apolparm->setgrid) {
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]);
1865  }
1866  spacs[i] = apolparm->grid[i];
1867  }
1868  }
1869 
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;
1876 
1877  for (x=xmin; x<=xmax; x=x+spacs[0]) {
1878  if ( VABS(x - xmin) < VSMALL) {
1879  wx = 0.5;
1880  } else if ( VABS(x - xmax) < VSMALL) {
1881  wx = 0.5;
1882  } else {
1883  wx = 1.0;
1884  }
1885  vec[0] = x;
1886  for (y=ymin; y<=ymax; y=y+spacs[1]) {
1887  if ( VABS(y - ymin) < VSMALL) {
1888  wy = 0.5;
1889  } else if ( VABS(y - ymax) < VSMALL) {
1890  wy = 0.5;
1891  } else {
1892  wy = 1.0;
1893  }
1894  vec[1] = y;
1895  for (z=zmin; z<=zmax; z=z+spacs[2]) {
1896  if ( VABS(z - zmin) < VSMALL) {
1897  wz = 0.5;
1898  } else if ( VABS(z - zmax) < VSMALL) {
1899  wz = 0.5;
1900  } else {
1901  wz = 1.0;
1902  }
1903  vec[2] = z;
1904 
1905  w = wx*wy*wz;
1906 
1907  chi = Vacc_ivdwAcc(thee, vec, srad);
1908 
1909  if (chi != 0.0) {
1910 
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);
1915 
1916  if (r <= 14 && r >= sigma){
1917 
1918  fo = 12.0*chi*rho*epsilon*(sigma6/VPOW(r,7)-sigma12/VPOW(r,13));
1919 
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;
1923 
1924  }else {
1925  for (si=0; si < 3; si++) fpt[si] = 0.0;
1926  }
1927  }else {
1928  for (si=0; si < 3; si++) fpt[si] = 0.0;
1929  }
1930 
1931  for(i=0;i<3;i++){
1932  force[i] += (w*fpt[i]);
1933  }
1934 
1935  } /* z loop */
1936  } /* y loop */
1937  } /* x loop */
1938 
1939  w = spacs[0]*spacs[1]*spacs[2];
1940  for(i=0;i<3;i++) force[i] *= w;
1941 
1942  return VRC_SUCCESS;
1943 }
1944 
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
Definition: vacc.c:132
double position[3]
Definition: vatom.h:86
double * xpts
Definition: vacc.h:86
double bconc
Definition: apolparm.h:139
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:1592
VPUBLIC void Vacc_totalAtomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Testing purposes only.
Definition: vacc.c:1401
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
Definition: vacc.c:780
double upper_corner[VAPBS_DIM]
Definition: vclist.h:127
Oracle for solvent- and ion-accessibility around a biomolecule.
Definition: vacc.h:108
VPRIVATE int Vacc_storeParms(Vacc *thee, Valist *alist, Vclist *clist, double surf_density)
Definition: vacc.c:148
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
Definition: vacc.c:245
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
VaccSurf ** surf
Definition: vacc.h:117
VPUBLIC double Vclist_maxRadius(Vclist *thee)
Get the max probe radius value (in A) the cell list was constructed with.
Definition: vclist.c:68
Contains declarations for class Vacc.
char * bpts
Definition: vacc.h:89
VPUBLIC unsigned long int Vacc_memChk(Vacc *thee)
Get number of bytes in this object and its members.
Definition: vacc.c:63
VaccSurf * refSphere
Definition: vacc.h:116
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
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...
Definition: vacc.c:377
double srad
Definition: apolparm.h:154
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition: valist.c:105
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:1018
VPUBLIC void VaccSurf_dtor2(VaccSurf *thee)
Destroy the surface object.
Definition: vacc.c:858
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
Definition: vacc.c:1733
VPUBLIC VaccSurf * VaccSurf_refSphere(Vmem *mem, int npts)
Set up an array of points for a reference sphere of unit radius.
Definition: vacc.c:938
Vmem * mem
Definition: vacc.h:85
int * atomFlags
Definition: vacc.h:113
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 Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
Definition: vacc.c:774
double probe_radius
Definition: vacc.h:93
int setgrid
Definition: apolparm.h:134
VPUBLIC int Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
Definition: vacc.c:1768
VPUBLIC void Vacc_dtor2(Vacc *thee)
FORTRAN stub to destroy object.
Definition: vacc.c:255
int natoms
Definition: vclist.h:103
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:1111
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
VEXTERNC double Vacc_vdwAcc(Vacc *thee, double center[VAPBS_DIM])
Report van der Waals accessibility.
Atom cell list cell.
Definition: vclist.h:101
double surf_density
Definition: vacc.h:122
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
double radius
Definition: vatom.h:87
double * zpts
Definition: vacc.h:88
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
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
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition: vhal.h:563
VPUBLIC double Vatom_getAtomID(Vatom *thee)
Get atom ID.
Definition: vatom.c:84
Vclist * clist
Definition: vacc.h:112
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
Definition: vacc.c:1515
VPUBLIC void Vacc_totalAtomdSAV(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA, Vclist *clist)
Total solvent accessible volume.
Definition: vacc.c:1460
double epsilon
Definition: vatom.h:91
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition: valist.c:115
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
double grid[3]
Definition: apolparm.h:133
int setwat
Definition: apolparm.h:180
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:818
#define VAPBS_DIM
Our dimension.
Definition: vhal.h:403
VPUBLIC VaccSurf * Vacc_atomSASPoints(Vacc *thee, double radius, Vatom *atom)
Get the set of points for this atom&#39;s solvent-accessible surface.
Definition: vacc.c:994
double watepsilon
Definition: apolparm.h:174
Surface object list of per-atom surface points.
Definition: vacc.h:84
double area
Definition: vacc.h:91
Vatom ** atoms
Definition: vclist.h:102
VPUBLIC void VaccSurf_dtor(VaccSurf **thee)
Destroy the surface object and free its memory.
Definition: vacc.c:844
VPRIVATE int Vacc_allocate(Vacc *thee)
Definition: vacc.c:193
VPUBLIC double Vacc_molAcc(Vacc *thee, double center[VAPBS_DIM], double radius)
Report molecular accessibility.
Definition: vacc.c:608
#define MAX_SPHERE_PTS
Maximum number of points on a sphere.
Definition: vhal.h:415
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition: vatom.c:63
VPUBLIC int Vacc_ctor2(Vacc *thee, Valist *alist, Vclist *clist, double surf_density)
FORTRAN stub to construct the accessibility object.
Definition: vacc.c:213
Contains public data members for Vatom class/module.
Definition: vatom.h:84
Container class for list of atom objects.
Definition: valist.h:78
double srad
Definition: pbeparm.h:152
double watsigma
Definition: apolparm.h:173
double lower_corner[VAPBS_DIM]
Definition: vclist.h:126
#define Vunit_pi
Pi.
Definition: vunit.h:104
Atom cell list.
Definition: vclist.h:117
int npts
Definition: vacc.h:92
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
Definition: vatom.c:105
int id
Definition: vatom.h:93
Vmem * mem
Definition: vacc.h:110
Valist * alist
Definition: vacc.h:111
double wcaEnergy
Definition: apolparm.h:177
VPUBLIC double Vacc_splineAcc(Vacc *thee, double center[VAPBS_DIM], double win, double infrad)
Report spline-based accessibility.
Definition: vacc.c:528
Parameter structure for APOL-specific variables from input files.
Definition: apolparm.h:129
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:873
double * ypts
Definition: vacc.h:87