APBS  1.5
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  thee = (VaccSurf*)calloc(1,sizeof(Vacc));
808  VASSERT( VaccSurf_ctor2(thee, mem, probe_radius, nsphere) );
809 
810  return thee;
811 }
812 
813 VPUBLIC 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 
839 VPUBLIC 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 
853 VPUBLIC 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 
868 VPUBLIC 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 
926 VPUBLIC 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 
982 VPUBLIC 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 
1006 VPUBLIC 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 
1099 VPUBLIC 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)
1200 VPUBLIC 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 */
1252 VPRIVATE 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 
1277 VPRIVATE 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)
1320 VPUBLIC 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 */
1389 VPUBLIC 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 */
1448 VPUBLIC 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 
1503 VPUBLIC 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 
1580 int 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 
1721 VPUBLIC 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 
1756 VPUBLIC 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 }
1932 
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:1580
VPUBLIC void Vacc_totalAtomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Testing purposes only.
Definition: vacc.c:1389
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:1006
VPUBLIC void VaccSurf_dtor2(VaccSurf *thee)
Destroy the surface object.
Definition: vacc.c:853
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
Definition: vacc.c:1721
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
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:1756
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:1099
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:556
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:1503
VPUBLIC void Vacc_totalAtomdSAV(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA, Vclist *clist)
Total solvent accessible volume.
Definition: vacc.c:1448
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:813
#define VAPBS_DIM
Our dimension.
Definition: vhal.h:402
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
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:839
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
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 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:868
double * ypts
Definition: vacc.h:87