APBS 3.0.0
Loading...
Searching...
No Matches
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
76# define FMM_DIST_TOL VSMALL
82# define FMM_IFLAG 2
86# define FMM_ORDER 4
90# define FMM_THETA 0.5
94# define FMM_MAXPARNODE 150
98# define FMM_SHRINK 1
102# define FMM_MINLEVEL 50000
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 */
117VPRIVATE 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 */
126VPRIVATE 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 */
135VPRIVATE 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
143
144 VASSERT(thee != VNULL);
145 return thee->alist;
146
147}
148
149VPUBLIC 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
156VPUBLIC Vgreen* Vgreen_ctor(Valist *alist) {
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
167VPUBLIC 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
192VPUBLIC 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
200VPUBLIC void Vgreen_dtor2(Vgreen *thee) {
201
202#ifdef HAVE_TREE
203 treecleanup(thee);
204#endif
205 Vmem_dtor(&(thee->vmem));
206
207}
208
209VPUBLIC 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
216VPUBLIC 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
224VPUBLIC 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
258VPUBLIC 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
310VPUBLIC 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
362VPUBLIC 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
425VPRIVATE 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
505VPRIVATE 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
525VPRIVATE 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}
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition valist.c:115
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition valist.c:105
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition vatom.c:63
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
Definition vatom.c:119
VPUBLIC void Vgreen_dtor(Vgreen **thee)
Destruct the Green's function oracle.
Definition vgreen.c:192
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's Law Green's function (solution to Laplace's equation) integrated over t...
Definition vgreen.c:310
VPUBLIC Valist * Vgreen_getValist(Vgreen *thee)
Get the atom list associated with this Green's function object.
Definition vgreen.c:142
VPUBLIC void Vgreen_dtor2(Vgreen *thee)
FORTRAN stub to destruct the Green's function oracle.
Definition vgreen.c:200
VPUBLIC int Vgreen_helmholtz(Vgreen *thee, int npos, double *x, double *y, double *z, double *val, double kappa)
Get the Green's function for Helmholtz's equation integrated over the atomic point charges.
Definition vgreen.c:209
VPUBLIC int Vgreen_ctor2(Vgreen *thee, Valist *alist)
FORTRAN stub to construct the Green's function oracle.
Definition vgreen.c:167
VPUBLIC Vgreen * Vgreen_ctor(Valist *alist)
Construct the Green's function oracle.
Definition vgreen.c:156
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's Law Green's function (solution to Laplace's equation) integrated over t...
Definition vgreen.c:362
VPUBLIC int Vgreen_coulomb(Vgreen *thee, int npos, double *x, double *y, double *z, double *val)
Get the Coulomb's Law Green's function (solution to Laplace's equation) integrated over the atomic po...
Definition vgreen.c:258
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 int Vgreen_coulomb_direct(Vgreen *thee, int npos, double *x, double *y, double *z, double *val)
Get the Coulomb's Law Green's function (solution to Laplace's equation) integrated over the atomic po...
Definition vgreen.c:224
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's function for Helmholtz's equation integrated over the atomic point charge...
Definition vgreen.c:216
#define Vunit_ec
Charge of an electron in C.
Definition vunit.h:92
#define Vunit_pi
Pi.
Definition vunit.h:104
#define Vunit_eps0
Vacuum permittivity.
Definition vunit.h:108
Container class for list of atom objects.
Definition valist.h:78
Contains public data members for Vatom class/module.
Definition vatom.h:84
Contains public data members for Vgreen class/module.
Definition vgreen.h:82
double * zp
Definition vgreen.h:90
double * xp
Definition vgreen.h:86
int np
Definition vgreen.h:94
double * yp
Definition vgreen.h:88
double * qp
Definition vgreen.h:92
Vmem * vmem
Definition vgreen.h:85
Valist * alist
Definition vgreen.h:84
Contains declarations for class Vgreen.