1 /**
2 * @file vgreen.c
3 * @ingroup Vgreen
4 * @author Nathan Baker
5 * @brief Class Vgreen methods
6 * @version $Id$
7 * @attention
8 * @verbatim
9 *
10 * APBS -- Adaptive Poisson-Boltzmann Solver
11 *
12 * Nathan A. Baker (nathan.baker@pnnl.gov)
13 * Pacific Northwest National Laboratory
14 *
15 * Additional contributing authors listed in the code documentation.
16 *
17 * Copyright (c) 2010-2014 Battelle Memorial Institute. Developed at the
18 * Pacific Northwest National Laboratory, operated by Battelle Memorial
19 * Institute, Pacific Northwest Division for the U.S. Department of Energy.
20 *
21 * Portions Copyright (c) 2002-2010, Washington University in St. Louis.
22 * Portions Copyright (c) 2002-2010, Nathan A. Baker.
23 * Portions Copyright (c) 1999-2002, The Regents of the University of
24 * California.
25 * Portions Copyright (c) 1995, Michael Holst.
26 * All rights reserved.
27 *
28 * Redistribution and use in source and binary forms, with or without
29 * modification, are permitted provided that the following conditions are met:
30 *
31 * Redistributions of source code must retain the above copyright notice, this
32 * list of conditions and the following disclaimer.
33 *
34 * Redistributions in binary form must reproduce the above copyright notice,
35 * this list of conditions and the following disclaimer in the documentation
36 * and/or other materials provided with the distribution.
37 *
38 * Neither the name of the developer nor the names of its contributors may be
39 * used to endorse or promote products derived from this software without
40 * specific prior written permission.
41 *
42 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
43 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
44 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
45 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
46 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
47 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
48 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
49 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
50 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
51 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
52 * THE POSSIBILITY OF SUCH DAMAGE.
53 *
54 * @endverbatim
55 */
56
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 /**
74 * @brief Lower distance cutoff for electrostatic interactions
75 * @ingroup Vgreen */
76 # define FMM_DIST_TOL VSMALL
77 /**
78 * @brief Flag for energy and force evaluation:
79 * \li 1 => evaluate energy only
80 * \li 2 => evaluate energy and force
81 * @ingroup Vgreen */
82 # define FMM_IFLAG 2
83 /**
84 * @brief Order of multipole expansion
85 * @ingroup Vgreen */
86 # define FMM_ORDER 4
87 /**
88 * @brief Multipole acceptance criterion
89 * @ingroup Vgreen */
90 # define FMM_THETA 0.5
91 /**
92 * @brief Maximum number of particles per node
93 * @ingroup Vgreen */
94 # define FMM_MAXPARNODE 150
95 /**
96 * @brief Switch for oct-tree construction
97 * @ingroup Vgreen */
98 # define FMM_SHRINK 1
99 /**
100 * @brief Minimum treecode level
101 * @ingroup Vgreen */
102 # define FMM_MINLEVEL 50000
103 /**
104 * @brief Maximum treecode level
105 * @ingroup Vgreen */
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
Vgreen_getValist(Vgreen * thee)142 VPUBLIC Valist* Vgreen_getValist(Vgreen *thee) {
143
144 VASSERT(thee != VNULL);
145 return thee->alist;
146
147 }
148
Vgreen_memChk(Vgreen * thee)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
Vgreen_ctor(Valist * alist)156 VPUBLIC 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
Vgreen_ctor2(Vgreen * thee,Valist * alist)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
Vgreen_dtor(Vgreen ** thee)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
Vgreen_dtor2(Vgreen * thee)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
Vgreen_helmholtz(Vgreen * thee,int npos,double * x,double * y,double * z,double * val,double kappa)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
Vgreen_helmholtzD(Vgreen * thee,int npos,double * x,double * y,double * z,double * gradx,double * grady,double * gradz,double kappa)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
Vgreen_coulomb_direct(Vgreen * thee,int npos,double * x,double * y,double * z,double * val)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
Vgreen_coulomb(Vgreen * thee,int npos,double * x,double * y,double * z,double * val)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
Vgreen_coulombD_direct(Vgreen * thee,int npos,double * x,double * y,double * z,double * pot,double * gradx,double * grady,double * gradz)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
Vgreen_coulombD(Vgreen * thee,int npos,double * x,double * y,double * z,double * pot,double * gradx,double * grady,double * gradz)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
treesetup(Vgreen * thee)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
treecleanup(Vgreen * thee)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
treecalc(Vgreen * thee,double * xtar,double * ytar,double * ztar,double * qtar,int numtars,double * tpengtar,double * x,double * y,double * z,double * q,int numpars,double * fx,double * fy,double * fz,int iflag,int farrdim,int arrdim)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
554