APBS  1.4.1
vgreen.c
Go to the documentation of this file.
1 
57 #include "vgreen.h"
58 
59 /* Define wrappers for F77 treecode routines */
60 #ifdef HAVE_TREE
61 # define F77TREEPEFORCE VF77_MANGLE(treepeforce, TREEPEFORCE)
62 # define F77DIRECT_ENG_FORCE VF77_MANGLE(direct_eng_force, DIRECT_ENG_FORCE)
63 # define F77CLEANUP VF77_MANGLE(mycleanup, MYCLEANUP)
64 # define F77TREE_COMPP VF77_MANGLE(mytree_compp, MYTREE_COMP)
65 # define F77TREE_COMPFP VF77_MANGLE(mytree_compfp, MYTREE_COMPFP)
66 # define F77CREATE_TREE VF77_MANGLE(mycreate_tree, MYCREATE_TREE)
67 # define F77INITLEVELS VF77_MANGLE(myinitlevels, MYINITLEVELS)
68 # define F77SETUP VF77_MANGLE(mysetup, MYSETUP)
69 #endif /* ifdef HAVE_TREE */
70 
71 /* Some constants associated with the tree code */
72 #ifdef HAVE_TREE
73 
76 # define FMM_DIST_TOL VSMALL
77 
82 # define FMM_IFLAG 2
83 
86 # define FMM_ORDER 4
87 
90 # define FMM_THETA 0.5
91 
94 # define FMM_MAXPARNODE 150
95 
98 # define FMM_SHRINK 1
99 
102 # define FMM_MINLEVEL 50000
103 
106 # define FMM_MAXLEVEL 0
107 #endif /* ifdef HAVE_TREE */
108 
109 
110 /*
111  * @brief Setup treecode internal structures
112  * @ingroup Vgreen
113  * @author Nathan Baker
114  * @param thee Vgreen object
115  * @return 1 if successful, 0 otherwise
116  */
117 VPRIVATE int treesetup(Vgreen *thee);
118 
119 /*
120  * @brief Clean up treecode internal structures
121  * @ingroup Vgreen
122  * @author Nathan Baker
123  * @param thee Vgreen object
124  * @return 1 if successful, 0 otherwise
125  */
126 VPRIVATE int treecleanup(Vgreen *thee);
127 
128 /*
129  * @brief Calculate forces or potential
130  * @ingroup Vgreen
131  * @author Nathan Baker
132  * @param thee Vgreen object
133  * @return 1 if successful, 0 otherwise
134  */
135 VPRIVATE int treecalc(Vgreen *thee, double *xtar, double *ytar, double *ztar,
136  double *qtar, int numtars, double *tpengtar, double *x, double *y,
137  double *z, double *q, int numpars, double *fx, double *fy, double *fz,
138  int iflag, int farrdim, int arrdim);
139 
140 #if !defined(VINLINE_VGREEN)
141 
142 VPUBLIC Valist* Vgreen_getValist(Vgreen *thee) {
143 
144  VASSERT(thee != VNULL);
145  return thee->alist;
146 
147 }
148 
149 VPUBLIC unsigned long int Vgreen_memChk(Vgreen *thee) {
150  if (thee == VNULL) return 0;
151  return Vmem_bytes(thee->vmem);
152 }
153 
154 #endif /* if !defined(VINLINE_VGREEN) */
155 
157 
158  /* Set up the structure */
159  Vgreen *thee = VNULL;
160  thee = (Vgreen *)Vmem_malloc(VNULL, 1, sizeof(Vgreen) );
161  VASSERT( thee != VNULL);
162  VASSERT( Vgreen_ctor2(thee, alist));
163 
164  return thee;
165 }
166 
167 VPUBLIC int Vgreen_ctor2(Vgreen *thee, Valist *alist) {
168 
169  VASSERT( thee != VNULL );
170 
171  /* Memory management object */
172  thee->vmem = Vmem_ctor("APBS:VGREEN");
173 
174  /* Set up the atom list and grid manager */
175  if (alist == VNULL) {
176  Vnm_print(2,"Vgreen_ctor2: got null pointer to Valist object!\n");
177  }
178 
179  thee->alist = alist;
180 
181  /* Setup FMM tree (if applicable) */
182 #ifdef HAVE_TREE
183  if (!treesetup(thee)) {
184  Vnm_print(2, "Vgreen_ctor2: Error setting up FMM tree!\n");
185  return 0;
186  }
187 #endif /* ifdef HAVE_TREE */
188 
189  return 1;
190 }
191 
192 VPUBLIC void Vgreen_dtor(Vgreen **thee) {
193  if ((*thee) != VNULL) {
194  Vgreen_dtor2(*thee);
195  Vmem_free(VNULL, 1, sizeof(Vgreen), (void **)thee);
196  (*thee) = VNULL;
197  }
198 }
199 
200 VPUBLIC void Vgreen_dtor2(Vgreen *thee) {
201 
202 #ifdef HAVE_TREE
203  treecleanup(thee);
204 #endif
205  Vmem_dtor(&(thee->vmem));
206 
207 }
208 
209 VPUBLIC int Vgreen_helmholtz(Vgreen *thee, int npos, double *x, double *y,
210  double *z, double *val, double kappa) {
211 
212  Vnm_print(2, "Error -- Vgreen_helmholtz not implemented yet!\n");
213  return 0;
214 }
215 
216 VPUBLIC int Vgreen_helmholtzD(Vgreen *thee, int npos, double *x, double *y,
217  double *z, double *gradx, double *grady, double *gradz, double kappa) {
218 
219  Vnm_print(2, "Error -- Vgreen_helmholtzD not implemented yet!\n");
220  return 0;
221 
222 }
223 
224 VPUBLIC int Vgreen_coulomb_direct(Vgreen *thee, int npos, double *x,
225  double *y, double *z, double *val) {
226 
227  Vatom *atom;
228  double *apos, charge, dist, dx, dy, dz, scale;
229  double *q, qtemp, fx, fy, fz;
230  int iatom, ipos;
231 
232  if (thee == VNULL) {
233  Vnm_print(2, "Vgreen_coulomb: Got NULL thee!\n");
234  return 0;
235  }
236 
237  for (ipos=0; ipos<npos; ipos++) val[ipos] = 0.0;
238 
239  for (iatom=0; iatom<Valist_getNumberAtoms(thee->alist); iatom++) {
240  atom = Valist_getAtom(thee->alist, iatom);
241  apos = Vatom_getPosition(atom);
242  charge = Vatom_getCharge(atom);
243  for (ipos=0; ipos<npos; ipos++) {
244  dx = apos[0] - x[ipos];
245  dy = apos[1] - y[ipos];
246  dz = apos[2] - z[ipos];
247  dist = VSQRT(VSQR(dx) + VSQR(dy) + VSQR(dz));
248  if (dist > VSMALL) val[ipos] += (charge/dist);
249  }
250  }
251 
252  scale = Vunit_ec/(4*Vunit_pi*Vunit_eps0*1.0e-10);
253  for (ipos=0; ipos<npos; ipos++) val[ipos] = val[ipos]*scale;
254 
255  return 1;
256 }
257 
258 VPUBLIC int Vgreen_coulomb(Vgreen *thee, int npos, double *x, double *y,
259  double *z, double *val) {
260 
261  Vatom *atom;
262  double *apos, charge, dist, dx, dy, dz, scale;
263  double *q, qtemp, fx, fy, fz;
264  int iatom, ipos;
265 
266  if (thee == VNULL) {
267  Vnm_print(2, "Vgreen_coulomb: Got NULL thee!\n");
268  return 0;
269  }
270 
271  for (ipos=0; ipos<npos; ipos++) val[ipos] = 0.0;
272 
273 #ifdef HAVE_TREE
274 
275  /* Allocate charge array (if necessary) */
276  if (Valist_getNumberAtoms(thee->alist) > 1) {
277  if (npos > 1) {
278  q = VNULL;
279  q = Vmem_malloc(thee->vmem, npos, sizeof(double));
280  if (q == VNULL) {
281  Vnm_print(2, "Vgreen_coulomb: Error allocating charge array!\n");
282  return 0;
283  }
284  } else {
285  q = &(qtemp);
286  }
287  for (ipos=0; ipos<npos; ipos++) q[ipos] = 1.0;
288 
289  /* Calculate */
290  treecalc(thee, x, y, z, q, npos, val, thee->xp, thee->yp, thee->zp,
291  thee->qp, thee->np, &fx, &fy, &fz, 1, 1, thee->np);
292  } else return Vgreen_coulomb_direct(thee, npos, x, y, z, val);
293 
294  /* De-allocate charge array (if necessary) */
295  if (npos > 1) Vmem_free(thee->vmem, npos, sizeof(double), (void **)&q);
296 
297  scale = Vunit_ec/(4*Vunit_pi*Vunit_eps0*1.0e-10);
298  for (ipos=0; ipos<npos; ipos++) val[ipos] = val[ipos]*scale;
299 
300  return 1;
301 
302 #else /* ifdef HAVE_TREE */
303 
304  return Vgreen_coulomb_direct(thee, npos, x, y, z, val);
305 
306 #endif
307 
308 }
309 
310 VPUBLIC int Vgreen_coulombD_direct(Vgreen *thee, int npos,
311  double *x, double *y, double *z, double *pot, double *gradx,
312  double *grady, double *gradz) {
313 
314  Vatom *atom;
315  double *apos, charge, dist, dist2, idist3, dy, dz, dx, scale;
316  double *q, qtemp;
317  int iatom, ipos;
318 
319  if (thee == VNULL) {
320  Vnm_print(2, "Vgreen_coulombD: Got VNULL thee!\n");
321  return 0;
322  }
323 
324  for (ipos=0; ipos<npos; ipos++) {
325  pot[ipos] = 0.0;
326  gradx[ipos] = 0.0;
327  grady[ipos] = 0.0;
328  gradz[ipos] = 0.0;
329  }
330 
331  for (iatom=0; iatom<Valist_getNumberAtoms(thee->alist); iatom++) {
332  atom = Valist_getAtom(thee->alist, iatom);
333  apos = Vatom_getPosition(atom);
334  charge = Vatom_getCharge(atom);
335  for (ipos=0; ipos<npos; ipos++) {
336  dx = apos[0] - x[ipos];
337  dy = apos[1] - y[ipos];
338  dz = apos[2] - z[ipos];
339  dist2 = VSQR(dx) + VSQR(dy) + VSQR(dz);
340  dist = VSQRT(dist2);
341  if (dist > VSMALL) {
342  idist3 = 1.0/(dist*dist2);
343  gradx[ipos] -= (charge*dx*idist3);
344  grady[ipos] -= (charge*dy*idist3);
345  gradz[ipos] -= (charge*dz*idist3);
346  pot[ipos] += (charge/dist);
347  }
348  }
349  }
350 
351  scale = Vunit_ec/(4*VPI*Vunit_eps0*(1.0e-10));
352  for (ipos=0; ipos<npos; ipos++) {
353  gradx[ipos] = gradx[ipos]*scale;
354  grady[ipos] = grady[ipos]*scale;
355  gradz[ipos] = gradz[ipos]*scale;
356  pot[ipos] = pot[ipos]*scale;
357  }
358 
359  return 1;
360 }
361 
362 VPUBLIC int Vgreen_coulombD(Vgreen *thee, int npos, double *x, double *y,
363  double *z, double *pot, double *gradx, double *grady, double *gradz) {
364 
365  Vatom *atom;
366  double *apos, charge, dist, dist2, idist3, dy, dz, dx, scale;
367  double *q, qtemp;
368  int iatom, ipos;
369 
370  if (thee == VNULL) {
371  Vnm_print(2, "Vgreen_coulombD: Got VNULL thee!\n");
372  return 0;
373  }
374 
375  for (ipos=0; ipos<npos; ipos++) {
376  pot[ipos] = 0.0;
377  gradx[ipos] = 0.0;
378  grady[ipos] = 0.0;
379  gradz[ipos] = 0.0;
380  }
381 
382 #ifdef HAVE_TREE
383 
384  if (Valist_getNumberAtoms(thee->alist) > 1) {
385  if (npos > 1) {
386  q = VNULL;
387  q = Vmem_malloc(thee->vmem, npos, sizeof(double));
388  if (q == VNULL) {
389  Vnm_print(2, "Vgreen_coulomb: Error allocating charge array!\n");
390  return 0;
391  }
392  } else {
393  q = &(qtemp);
394  }
395  for (ipos=0; ipos<npos; ipos++) q[ipos] = 1.0;
396 
397  /* Calculate */
398  treecalc(thee, x, y, z, q, npos, pot, thee->xp, thee->yp, thee->zp,
399  thee->qp, thee->np, gradx, grady, gradz, 2, npos, thee->np);
400 
401  /* De-allocate charge array (if necessary) */
402  if (npos > 1) Vmem_free(thee->vmem, npos, sizeof(double), (void **)&q);
403  } else return Vgreen_coulombD_direct(thee, npos, x, y, z, pot,
404  gradx, grady, gradz);
405 
406  scale = Vunit_ec/(4*VPI*Vunit_eps0*(1.0e-10));
407  for (ipos=0; ipos<npos; ipos++) {
408  gradx[ipos] = gradx[ipos]*scale;
409  grady[ipos] = grady[ipos]*scale;
410  gradz[ipos] = gradz[ipos]*scale;
411  pot[ipos] = pot[ipos]*scale;
412  }
413 
414  return 1;
415 
416 #else /* ifdef HAVE_TREE */
417 
418  return Vgreen_coulombD_direct(thee, npos, x, y, z, pot,
419  gradx, grady, gradz);
420 
421 #endif
422 
423 }
424 
425 VPRIVATE int treesetup(Vgreen *thee) {
426 
427 #ifdef HAVE_TREE
428 
429  double dist_tol = FMM_DIST_TOL;
430  int iflag = FMM_IFLAG;
431  double order = FMM_ORDER;
432  int theta = FMM_THETA;
433  int shrink = FMM_SHRINK;
434  int maxparnode = FMM_MAXPARNODE;
435  int minlevel = FMM_MINLEVEL;
436  int maxlevel = FMM_MAXLEVEL;
437  int level = 0;
438  int one = 1;
439  Vatom *atom;
440  double xyzminmax[6], *pos;
441  int i;
442 
443  /* Set up particle arrays with atomic coordinates and charges */
444  Vnm_print(0, "treesetup: Initializing FMM particle arrays...\n");
445  thee->np = Valist_getNumberAtoms(thee->alist);
446  thee->xp = VNULL;
447  thee->xp = (double *)Vmem_malloc(thee->vmem, thee->np, sizeof(double));
448  if (thee->xp == VNULL) {
449  Vnm_print(2, "Vgreen_ctor2: Failed to allocate %d*sizeof(double)!\n",
450  thee->np);
451  return 0;
452  }
453  thee->yp = VNULL;
454  thee->yp = (double *)Vmem_malloc(thee->vmem, thee->np, sizeof(double));
455  if (thee->yp == VNULL) {
456  Vnm_print(2, "Vgreen_ctor2: Failed to allocate %d*sizeof(double)!\n",
457  thee->np);
458  return 0;
459  }
460  thee->zp = VNULL;
461  thee->zp = (double *)Vmem_malloc(thee->vmem, thee->np, sizeof(double));
462  if (thee->zp == VNULL) {
463  Vnm_print(2, "Vgreen_ctor2: Failed to allocate %d*sizeof(double)!\n",
464  thee->np);
465  return 0;
466  }
467  thee->qp = VNULL;
468  thee->qp = (double *)Vmem_malloc(thee->vmem, thee->np, sizeof(double));
469  if (thee->qp == VNULL) {
470  Vnm_print(2, "Vgreen_ctor2: Failed to allocate %d*sizeof(double)!\n",
471  thee->np);
472  return 0;
473  }
474  for (i=0; i<thee->np; i++) {
475  atom = Valist_getAtom(thee->alist, i);
476  pos = Vatom_getPosition(atom);
477  thee->xp[i] = pos[0];
478  thee->yp[i] = pos[1];
479  thee->zp[i] = pos[2];
480  thee->qp[i] = Vatom_getCharge(atom);
481  }
482 
483  Vnm_print(0, "treesetup: Setting things up...\n");
484  F77SETUP(thee->xp, thee->yp, thee->zp, &(thee->np), &order, &theta, &iflag,
485  &dist_tol, xyzminmax, &(thee->np));
486 
487 
488  Vnm_print(0, "treesetup: Initializing levels...\n");
489  F77INITLEVELS(&minlevel, &maxlevel);
490 
491  Vnm_print(0, "treesetup: Creating tree...\n");
492  F77CREATE_TREE(&one, &(thee->np), thee->xp, thee->yp, thee->zp, thee->qp,
493  &shrink, &maxparnode, xyzminmax, &level, &(thee->np));
494 
495  return 1;
496 
497 #else /* ifdef HAVE_TREE */
498 
499  Vnm_print(2, "treesetup: Error! APBS not linked with treecode!\n");
500  return 0;
501 
502 #endif /* ifdef HAVE_TREE */
503 }
504 
505 VPRIVATE int treecleanup(Vgreen *thee) {
506 
507 #ifdef HAVE_TREE
508 
509  Vmem_free(thee->vmem, thee->np, sizeof(double), (void **)&(thee->xp));
510  Vmem_free(thee->vmem, thee->np, sizeof(double), (void **)&(thee->yp));
511  Vmem_free(thee->vmem, thee->np, sizeof(double), (void **)&(thee->zp));
512  Vmem_free(thee->vmem, thee->np, sizeof(double), (void **)&(thee->qp));
513  F77CLEANUP();
514 
515  return 1;
516 
517 #else /* ifdef HAVE_TREE */
518 
519  Vnm_print(2, "treecleanup: Error! APBS not linked with treecode!\n");
520  return 0;
521 
522 #endif /* ifdef HAVE_TREE */
523 }
524 
525 VPRIVATE int treecalc(Vgreen *thee, double *xtar, double *ytar, double *ztar,
526  double *qtar, int numtars, double *tpengtar, double *x, double *y,
527  double *z, double *q, int numpars, double *fx, double *fy, double *fz,
528  int iflag, int farrdim, int arrdim) {
529 
530 #ifdef HAVE_TREE
531  int i, level, err, maxlevel, minlevel, one;
532  double xyzminmax[6];
533 
534 
535  if (iflag != 1) {
536  F77TREE_COMPFP(xtar, ytar, ztar, qtar, &numtars, tpengtar, x, y, z, q,
537  fx, fy, fz, &numpars, &farrdim, &arrdim);
538  } else {
539  F77TREE_COMPP(xtar, ytar, ztar, qtar, &numtars, tpengtar, &farrdim, x,
540  y, z, q, &numpars, &arrdim);
541  }
542 
543 
544  return 1;
545 
546 #else /* ifdef HAVE_TREE */
547 
548  Vnm_print(2, "treecalc: Error! APBS not linked with treecode!\n");
549  return 0;
550 
551 #endif /* ifdef HAVE_TREE */
552 }
553 
VPUBLIC unsigned long int Vgreen_memChk(Vgreen *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition: vgreen.c:149
VPUBLIC void Vgreen_dtor2(Vgreen *thee)
FORTRAN stub to destruct the Green&#39;s function oracle.
Definition: vgreen.c:200
Contains public data members for Vgreen class/module.
Definition: vgreen.h:82
Contains declarations for class Vgreen.
VPUBLIC int Vgreen_coulomb_direct(Vgreen *thee, int npos, double *x, double *y, double *z, double *val)
Get the Coulomb&#39;s Law Green&#39;s function (solution to Laplace&#39;s equation) integrated over the atomic po...
Definition: vgreen.c:224
double * zp
Definition: vgreen.h:90
VPUBLIC int Vgreen_helmholtzD(Vgreen *thee, int npos, double *x, double *y, double *z, double *gradx, double *grady, double *gradz, double kappa)
Get the gradient of Green&#39;s function for Helmholtz&#39;s equation integrated over the atomic point charge...
Definition: vgreen.c:216
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition: valist.c:105
VPUBLIC int Vgreen_coulomb(Vgreen *thee, int npos, double *x, double *y, double *z, double *val)
Get the Coulomb&#39;s Law Green&#39;s function (solution to Laplace&#39;s equation) integrated over the atomic po...
Definition: vgreen.c:258
VPUBLIC Valist * Vgreen_getValist(Vgreen *thee)
Get the atom list associated with this Green&#39;s function object.
Definition: vgreen.c:142
VPUBLIC Vgreen * Vgreen_ctor(Valist *alist)
Construct the Green&#39;s function oracle.
Definition: vgreen.c:156
Valist * alist
Definition: vgreen.h:84
VPUBLIC int Vgreen_ctor2(Vgreen *thee, Valist *alist)
FORTRAN stub to construct the Green&#39;s function oracle.
Definition: vgreen.c:167
VPUBLIC int Vgreen_coulombD_direct(Vgreen *thee, int npos, double *x, double *y, double *z, double *pot, double *gradx, double *grady, double *gradz)
Get gradient of the Coulomb&#39;s Law Green&#39;s function (solution to Laplace&#39;s equation) integrated over t...
Definition: vgreen.c:310
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition: valist.c:115
int np
Definition: vgreen.h:94
Valist * alist
Definition: vclist.h:120
VPUBLIC int Vgreen_coulombD(Vgreen *thee, int npos, double *x, double *y, double *z, double *pot, double *gradx, double *grady, double *gradz)
Get gradient of the Coulomb&#39;s Law Green&#39;s function (solution to Laplace&#39;s equation) integrated over t...
Definition: vgreen.c:362
Vmem * vmem
Definition: vgreen.h:85
VPUBLIC int Vgreen_helmholtz(Vgreen *thee, int npos, double *x, double *y, double *z, double *val, double kappa)
Get the Green&#39;s function for Helmholtz&#39;s equation integrated over the atomic point charges...
Definition: vgreen.c:209
double * qp
Definition: vgreen.h:92
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition: vatom.c:63
Contains public data members for Vatom class/module.
Definition: vatom.h:84
VPUBLIC void Vgreen_dtor(Vgreen **thee)
Destruct the Green&#39;s function oracle.
Definition: vgreen.c:192
Container class for list of atom objects.
Definition: valist.h:78
#define Vunit_eps0
Vacuum permittivity.
Definition: vunit.h:108
#define Vunit_pi
Pi.
Definition: vunit.h:104
double * xp
Definition: vgreen.h:86
double * yp
Definition: vgreen.h:88
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
Definition: vatom.c:119
#define Vunit_ec
Charge of an electron in C.
Definition: vunit.h:92