APBS 3.0.0
Loading...
Searching...
No Matches
vpbe.c
Go to the documentation of this file.
1
57#include "vpbe.h"
58
59/* ///////////////////////////////////////////////////////////////////////////
60// Class Vpbe: Private method declaration
62#define MAX_SPLINE_WINDOW 0.5
63
64/* ///////////////////////////////////////////////////////////////////////////
65// Class Vpbe: Inlineable methods
67#if !defined(VINLINE_VPBE)
68
69VPUBLIC Valist* Vpbe_getValist(Vpbe *thee) {
70
71 VASSERT(thee != VNULL);
72 return thee->alist;
73
74}
75
76VPUBLIC Vacc* Vpbe_getVacc(Vpbe *thee) {
77
78 VASSERT(thee != VNULL);
79 VASSERT(thee->paramFlag);
80 return thee->acc;
81
82}
83
84VPUBLIC double Vpbe_getBulkIonicStrength(Vpbe *thee) {
85
86 VASSERT(thee != VNULL);
87 VASSERT(thee->paramFlag);
88 return thee->bulkIonicStrength;
89}
90
91VPUBLIC double Vpbe_getTemperature(Vpbe *thee) {
92
93 VASSERT(thee != VNULL);
94 VASSERT(thee->paramFlag);
95 return thee->T;
96
97}
98
99VPUBLIC double Vpbe_getSoluteDiel(Vpbe *thee) {
100
101 VASSERT(thee != VNULL);
102 VASSERT(thee->paramFlag);
103 return thee->soluteDiel;
104
105}
106
107VPUBLIC double* Vpbe_getSoluteCenter(Vpbe *thee) {
108
109 VASSERT(thee != VNULL);
110 return thee->soluteCenter;
111}
112
113VPUBLIC double Vpbe_getSolventDiel(Vpbe *thee) {
114
115 VASSERT(thee != VNULL);
116 VASSERT(thee->paramFlag);
117 return thee->solventDiel;
118}
119
120VPUBLIC double Vpbe_getSolventRadius(Vpbe *thee) {
121
122 VASSERT(thee != VNULL);
123 VASSERT(thee->paramFlag);
124 return thee->solventRadius;
125}
126
127VPUBLIC double Vpbe_getMaxIonRadius(Vpbe *thee) {
128
129 VASSERT(thee != VNULL);
130 VASSERT(thee->paramFlag);
131 return thee->maxIonRadius;
132}
133
134VPUBLIC double Vpbe_getXkappa(Vpbe *thee) {
135
136 VASSERT(thee != VNULL);
137 VASSERT(thee->paramFlag);
138 return thee->xkappa;
139}
140
141VPUBLIC double Vpbe_getDeblen(Vpbe *thee) {
142
143 VASSERT(thee != VNULL);
144 VASSERT(thee->paramFlag);
145 return thee->deblen;
146}
147
148VPUBLIC double Vpbe_getZkappa2(Vpbe *thee) {
149
150 VASSERT(thee != VNULL);
151 VASSERT(thee->paramFlag);
152 return thee->zkappa2;
153}
154
155VPUBLIC double Vpbe_getZmagic(Vpbe *thee) {
156
157 VASSERT(thee != VNULL);
158 VASSERT(thee->paramFlag);
159 return thee->zmagic;
160}
161
162VPUBLIC double Vpbe_getSoluteRadius(Vpbe *thee) {
163
164 VASSERT(thee != VNULL);
165 return thee->soluteRadius;
166}
167
168VPUBLIC double Vpbe_getSoluteXlen(Vpbe *thee) {
169
170 VASSERT(thee != VNULL);
171 return thee->soluteXlen;
172}
173
174VPUBLIC double Vpbe_getSoluteYlen(Vpbe *thee) {
175
176 VASSERT(thee != VNULL);
177 return thee->soluteYlen;
178}
179
180VPUBLIC double Vpbe_getSoluteZlen(Vpbe *thee) {
181
182 VASSERT(thee != VNULL);
183 return thee->soluteZlen;
184}
185
186VPUBLIC double Vpbe_getSoluteCharge(Vpbe *thee) {
187
188 VASSERT(thee != VNULL);
189 return thee->soluteCharge;
190}
191
192/* ///////////////////////////////////////////////////////////////////////////
193 // Routine: Vpbe_getzmem
194 // Purpose: This routine returns values stored in the structure thee.
195 // Author: Michael Grabe
197VPUBLIC double Vpbe_getzmem(Vpbe *thee) {
198
199 VASSERT(thee != VNULL);
200 VASSERT(thee->param2Flag);
201 return thee->z_mem;
202}
203
204/* ///////////////////////////////////////////////////////////////////////////
205 // Routine: Vpbe_getLmem
206 // Purpose: This routine returns values stored in the structure thee.
207 // Author: Michael Grabe
209VPUBLIC double Vpbe_getLmem(Vpbe *thee) {
210
211 VASSERT(thee != VNULL);
212 VASSERT(thee->param2Flag);
213 return thee->L;
214}
215
216/* ///////////////////////////////////////////////////////////////////////////
217 // Routine: Vpbe_getmembraneDiel
218 // Purpose: This routine returns values stored in the structure thee.
219 // Author: Michael Grabe
221VPUBLIC double Vpbe_getmembraneDiel(Vpbe *thee) {
222
223 VASSERT(thee != VNULL);
224 VASSERT(thee->param2Flag);
225 return thee->membraneDiel;
226}
227
228/* ///////////////////////////////////////////////////////////////////////////
229 // Routine: Vpbe_getmemv
230 // Purpose: This routine returns values stored in the structure thee.
231 // Author: Michael Grabe
233VPUBLIC double Vpbe_getmemv(Vpbe *thee) {
234
235 VASSERT(thee != VNULL);
236 VASSERT(thee->param2Flag);
237 return thee->V;
238}
239
240#endif /* if !defined(VINLINE_VPBE) */
241
242/* ///////////////////////////////////////////////////////////////////////////
243// Class Vpbe: Non-inlineable methods
245
246VPUBLIC Vpbe* Vpbe_ctor(Valist *alist, int ionNum, double *ionConc,
247 double *ionRadii, double *ionQ, double T,
248 double soluteDiel, double solventDiel,
249 double solventRadius, int focusFlag, double sdens,
250 double z_mem, double L, double membraneDiel, double V ) {
251
252 /* Set up the structure */
253 Vpbe *thee = VNULL;
254 thee = (Vpbe*)Vmem_malloc(VNULL, 1, sizeof(Vpbe) );
255 VASSERT( thee != VNULL);
256 VASSERT( Vpbe_ctor2(thee, alist, ionNum, ionConc, ionRadii, ionQ,
257 T, soluteDiel, solventDiel, solventRadius, focusFlag, sdens,
258 z_mem, L, membraneDiel, V) );
259
260 return thee;
261}
262
263
264VPUBLIC int Vpbe_ctor2(Vpbe *thee, Valist *alist, int ionNum,
265 double *ionConc, double *ionRadii,
266 double *ionQ, double T, double soluteDiel,
267 double solventDiel, double solventRadius, int focusFlag,
268 double sdens, double z_mem, double L, double membraneDiel,
269 double V) {
270
271 int i, iatom, inhash[3];
272 double atomRadius;
273 Vatom *atom;
274 double center[3] = {0.0, 0.0, 0.0};
275 double lower_corner[3] = {0.0, 0.0, 0.0};
276 double upper_corner[3] = {0.0, 0.0, 0.0};
277 double disp[3], dist, radius, charge, xmin, xmax, ymin, ymax, zmin, zmax;
278 double x, y, z, netCharge;
279 double nhash[3];
280 const double N_A = 6.022045000e+23;
281 const double e_c = 4.803242384e-10;
282 const double k_B = 1.380662000e-16;
283 const double pi = 4. * VATAN(1.);
284
285 /* Set up memory management object */
286 thee->vmem = Vmem_ctor("APBS::VPBE");
287
288 VASSERT(thee != VNULL);
289 if (alist == VNULL) {
290 Vnm_print(2, "Vpbe_ctor2: Got null pointer to Valist object!\n");
291 return 0;
292 }
293
294 /* **** STUFF THAT GETS DONE FOR EVERYONE **** */
295 /* Set pointers */
296 thee->alist = alist;
297 thee->paramFlag = 0;
298
299 /* Determine solute center */
300 center[0] = thee->alist->center[0];
301 center[1] = thee->alist->center[1];
302 center[2] = thee->alist->center[2];
303 thee->soluteCenter[0] = center[0];
304 thee->soluteCenter[1] = center[1];
305 thee->soluteCenter[2] = center[2];
306
307 /* Determine solute length and charge*/
308 radius = 0;
309 atom = Valist_getAtom(thee->alist, 0);
310 xmin = Vatom_getPosition(atom)[0];
311 xmax = Vatom_getPosition(atom)[0];
312 ymin = Vatom_getPosition(atom)[1];
313 ymax = Vatom_getPosition(atom)[1];
314 zmin = Vatom_getPosition(atom)[2];
315 zmax = Vatom_getPosition(atom)[2];
316 charge = 0;
317 for (iatom=0; iatom<Valist_getNumberAtoms(thee->alist); iatom++) {
318 atom = Valist_getAtom(thee->alist, iatom);
319 atomRadius = Vatom_getRadius(atom);
320 x = Vatom_getPosition(atom)[0];
321 y = Vatom_getPosition(atom)[1];
322 z = Vatom_getPosition(atom)[2];
323 if ((x+atomRadius) > xmax) xmax = x + atomRadius;
324 if ((x-atomRadius) < xmin) xmin = x - atomRadius;
325 if ((y+atomRadius) > ymax) ymax = y + atomRadius;
326 if ((y-atomRadius) < ymin) ymin = y - atomRadius;
327 if ((z+atomRadius) > zmax) zmax = z + atomRadius;
328 if ((z-atomRadius) < zmin) zmin = z - atomRadius;
329 disp[0] = (x - center[0]);
330 disp[1] = (y - center[1]);
331 disp[2] = (z - center[2]);
332 dist = (disp[0]*disp[0]) + (disp[1]*disp[1]) + (disp[2]*disp[2]);
333 dist = VSQRT(dist) + atomRadius;
334 if (dist > radius) radius = dist;
335 charge += Vatom_getCharge(Valist_getAtom(thee->alist, iatom));
336 }
337 thee->soluteRadius = radius;
338 Vnm_print(0, "Vpbe_ctor2: solute radius = %g\n", radius);
339 thee->soluteXlen = xmax - xmin;
340 thee->soluteYlen = ymax - ymin;
341 thee->soluteZlen = zmax - zmin;
342 Vnm_print(0, "Vpbe_ctor2: solute dimensions = %g x %g x %g\n",
343 thee->soluteXlen, thee->soluteYlen, thee->soluteZlen);
344 thee->soluteCharge = charge;
345 Vnm_print(0, "Vpbe_ctor2: solute charge = %g\n", charge);
346
347 /* Set parameters */
348 thee->numIon = ionNum;
349 if (thee->numIon >= MAXION) {
350 Vnm_print(2, "Vpbe_ctor2: Too many ion species (MAX = %d)!\n",
351 MAXION);
352 return 0;
353 }
354 thee->bulkIonicStrength = 0.0;
355 thee->maxIonRadius = 0.0;
356 netCharge = 0.0;
357 for (i=0; i<thee->numIon; i++) {
358 thee->ionConc[i] = ionConc[i];
359 thee->ionRadii[i] = ionRadii[i];
360 if (ionRadii[i] > thee->maxIonRadius) thee->maxIonRadius = ionRadii[i];
361 thee->ionQ[i] = ionQ[i];
362 thee->bulkIonicStrength += (0.5*ionConc[i]*VSQR(ionQ[i]));
363 netCharge += (ionConc[i]*ionQ[i]);
364 }
365#ifndef VAPBSQUIET
366 Vnm_print(1, " Vpbe_ctor: Using max ion radius (%g A) for exclusion \
367function\n", thee->maxIonRadius);
368#endif
369 if (VABS(netCharge) > VSMALL) {
370 Vnm_print(2, "Vpbe_ctor2: You have a counterion charge imbalance!\n");
371 Vnm_print(2, "Vpbe_ctor2: Net charge conc. = %g M\n", netCharge);
372 return 0;
373 }
374 thee->T = T;
375 thee->soluteDiel = soluteDiel;
376 thee->solventDiel = solventDiel;
377 thee->solventRadius = solventRadius;
378
379 /* Compute parameters:
380 *
381 * kappa^2 = (8 pi N_A e_c^2) I_s / (1000 eps_w k_B T)
382 * kappa = 0.325567 * I_s^{1/2} angstroms^{-1}
383 * deblen = 1 / kappa
384 * = 3.071564378 * I_s^{1/2} angstroms
385 * \bar{kappa}^2 = eps_w * kappa^2
386 * zmagic = (4 * pi * e_c^2) / (k_B T) (we scale the diagonal later)
387 * = 7046.528838
388 */
389 if (thee->T == 0.0) {
390 Vnm_print(2, "Vpbe_ctor2: You set the temperature to 0 K.\n");
391 Vnm_print(2, "Vpbe_ctor2: That violates the 3rd Law of Thermo.!");
392 return 0;
393 }
394 if (thee->bulkIonicStrength == 0.) {
395 thee->xkappa = 0.;
396 thee->deblen = 0.;
397 thee->zkappa2 = 0.;
398 } else {
399 thee->xkappa = VSQRT( thee->bulkIonicStrength * 1.0e-16 *
400 ((8.0 * pi * N_A * e_c*e_c) /
401 (1000.0 * thee->solventDiel * k_B * T))
402 );
403 thee->deblen = 1. / thee->xkappa;
404 thee->zkappa2 = thee->solventDiel * VSQR(thee->xkappa);
405 }
406 Vnm_print(0, "Vpbe_ctor2: bulk ionic strength = %g\n",
407 thee->bulkIonicStrength);
408 Vnm_print(0, "Vpbe_ctor2: xkappa = %g\n", thee->xkappa);
409 Vnm_print(0, "Vpbe_ctor2: Debye length = %g\n", thee->deblen);
410 Vnm_print(0, "Vpbe_ctor2: zkappa2 = %g\n", thee->zkappa2);
411 thee->zmagic = ((4.0 * pi * e_c*e_c) / (k_B * thee->T)) * 1.0e+8;
412 Vnm_print(0, "Vpbe_ctor2: zmagic = %g\n", thee->zmagic);
413
414 /* Compute accessibility objects:
415 * - Allow for extra room in the case of spline windowing
416 * - Place some limits on the size of the hash table in the case of very
417 * large molecules
418 */
419 if (thee->maxIonRadius > thee->solventRadius)
420 radius = thee->maxIonRadius + MAX_SPLINE_WINDOW;
421 else radius = thee->solventRadius + MAX_SPLINE_WINDOW;
422
423 nhash[0] = (thee->soluteXlen)/0.5;
424 nhash[1] = (thee->soluteYlen)/0.5;
425 nhash[2] = (thee->soluteZlen)/0.5;
426 for (i=0; i<3; i++) inhash[i] = (int)(nhash[i]);
427
428 for (i=0;i<3;i++){
429 if (inhash[i] < 3) inhash[i] = 3;
430 if (inhash[i] > MAX_HASH_DIM) inhash[i] = MAX_HASH_DIM;
431 }
432 Vnm_print(0, "Vpbe_ctor2: Constructing Vclist with %d x %d x %d table\n",
433 inhash[0], inhash[1], inhash[2]);
434
435 thee->clist = Vclist_ctor(thee->alist, radius, inhash,
436 CLIST_AUTO_DOMAIN, lower_corner, upper_corner);
437
438 VASSERT(thee->clist != VNULL);
439 thee->acc = Vacc_ctor(thee->alist, thee->clist, sdens);
440
441 VASSERT(thee->acc != VNULL);
442
443 /* SMPBE Added */
444 thee->smsize = 0.0;
445 thee->smvolume = 0.0;
446 thee->ipkey = 0;
447
448 thee->paramFlag = 1;
449
450 /*-----------------------------------------------------------*/
451 /* added by Michael Grabe */
452 /*-----------------------------------------------------------*/
453
454 thee->z_mem = z_mem;
455 thee->L = L;
456 thee->membraneDiel = membraneDiel;
457 thee->V = V;
458
459 if (V != 0.0) thee->param2Flag = 1;
460 else thee->param2Flag = 0;
461
462 /*-----------------------------------------------------------*/
463
464 return 1;
465}
466
467VPUBLIC void Vpbe_dtor(Vpbe **thee) {
468 if ((*thee) != VNULL) {
469 Vpbe_dtor2(*thee);
470 Vmem_free(VNULL, 1, sizeof(Vpbe), (void **)thee);
471 (*thee) = VNULL;
472 }
473}
474
475VPUBLIC void Vpbe_dtor2(Vpbe *thee) {
476 Vclist_dtor(&(thee->clist));
477 Vacc_dtor(&(thee->acc));
478 Vmem_dtor(&(thee->vmem));
479}
480
481VPUBLIC double Vpbe_getCoulombEnergy1(Vpbe *thee) {
482
483 int i, j, k, natoms;
484
485 double dist, *ipos, *jpos, icharge, jcharge;
486 double energy = 0.0;
487 double eps, T;
488 Vatom *iatom, *jatom;
489 Valist *alist;
490
491 VASSERT(thee != VNULL);
492 alist = Vpbe_getValist(thee);
493 VASSERT(alist != VNULL);
494 natoms = Valist_getNumberAtoms(alist);
495
496 /* Do the sum */
497 for (i=0; i<natoms; i++) {
498 iatom = Valist_getAtom(alist,i);
499 icharge = Vatom_getCharge(iatom);
500 ipos = Vatom_getPosition(iatom);
501 for (j=i+1; j<natoms; j++) {
502 jatom = Valist_getAtom(alist,j);
503 jcharge = Vatom_getCharge(jatom);
504 jpos = Vatom_getPosition(jatom);
505 dist = 0;
506 for (k=0; k<3; k++) dist += ((ipos[k]-jpos[k])*(ipos[k]-jpos[k]));
507 dist = VSQRT(dist);
508 energy = energy + icharge*jcharge/dist;
509 }
510 }
511
512 /* Convert the result to J */
513 T = Vpbe_getTemperature(thee);
514 eps = Vpbe_getSoluteDiel(thee);
515 energy = energy*Vunit_ec*Vunit_ec/(4*Vunit_pi*Vunit_eps0*eps*(1.0e-10));
516
517 /* Scale by Boltzmann energy */
518 energy = energy/(Vunit_kb*T);
519
520 return energy;
521}
522
523VPUBLIC unsigned long int Vpbe_memChk(Vpbe *thee) {
524
525 unsigned long int memUse = 0;
526
527 if (thee == VNULL) return 0;
528
529 memUse = memUse + sizeof(Vpbe);
530 memUse = memUse + (unsigned long int)Vacc_memChk(thee->acc);
531
532 return memUse;
533}
534
535VPUBLIC int Vpbe_getIons(Vpbe *thee, int *nion, double ionConc[MAXION],
536 double ionRadii[MAXION], double ionQ[MAXION]) {
537
538 int i;
539
540 VASSERT(thee != VNULL);
541
542 *nion = thee->numIon;
543 for (i=0; i<(*nion); i++) {
544 ionConc[i] = thee->ionConc[i];
545 ionRadii[i] = thee->ionRadii[i];
546 ionQ[i] = thee->ionQ[i];
547 }
548
549 return *nion;
550}
VPUBLIC unsigned long int Vacc_memChk(Vacc *thee)
Get number of bytes in this object and its members.
Definition vacc.c:63
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
Definition vacc.c:245
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
Definition vacc.c:132
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_getRadius(Vatom *thee)
Get atomic position.
Definition vatom.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 Vclist_dtor(Vclist **thee)
Destroy object.
Definition vclist.c:397
VPUBLIC Vclist * Vclist_ctor(Valist *alist, double max_radius, int npts[VAPBS_DIM], Vclist_DomainMode mode, double lower_corner[VAPBS_DIM], double upper_corner[VAPBS_DIM])
Construct the cell list object.
Definition vclist.c:75
@ CLIST_AUTO_DOMAIN
Definition vclist.h:83
#define MAXION
The maximum number of ion species that can be involved in a single PBE calculation.
Definition vhal.h:377
struct sVpbe Vpbe
Declaration of the Vpbe class as the Vpbe structure.
Definition vpbe.h:144
VPUBLIC void Vpbe_dtor2(Vpbe *thee)
FORTRAN stub object destructor.
Definition vpbe.c:475
VPUBLIC double Vpbe_getCoulombEnergy1(Vpbe *thee)
Calculate coulombic energy of set of charges.
Definition vpbe.c:481
VPUBLIC void Vpbe_dtor(Vpbe **thee)
Object destructor.
Definition vpbe.c:467
VPUBLIC double Vpbe_getSoluteDiel(Vpbe *thee)
Get solute dielectric constant.
Definition vpbe.c:99
VPUBLIC unsigned long int Vpbe_memChk(Vpbe *thee)
Return the memory used by this structure (and its contents) in bytes.
Definition vpbe.c:523
VPUBLIC int Vpbe_ctor2(Vpbe *thee, Valist *alist, int ionNum, double *ionConc, double *ionRadii, double *ionQ, double T, double soluteDiel, double solventDiel, double solventRadius, int focusFlag, double sdens, double z_mem, double L, double membraneDiel, double V)
FORTRAN stub to construct Vpbe objct.
Definition vpbe.c:264
VPUBLIC Valist * Vpbe_getValist(Vpbe *thee)
Get atom list.
Definition vpbe.c:69
VPUBLIC int Vpbe_getIons(Vpbe *thee, int *nion, double ionConc[MAXION], double ionRadii[MAXION], double ionQ[MAXION])
Get information about the counterion species present.
Definition vpbe.c:535
VPUBLIC double Vpbe_getTemperature(Vpbe *thee)
Get temperature.
Definition vpbe.c:91
#define Vunit_ec
Charge of an electron in C.
Definition vunit.h:92
#define Vunit_kb
Boltzmann constant.
Definition vunit.h:96
#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
double center[3]
Definition valist.h:81
Contains public data members for Vatom class/module.
Definition vatom.h:84
Contains public data members for Vpbe class/module.
Definition vpbe.h:84
double smsize
Definition vpbe.h:121
double soluteYlen
Definition vpbe.h:116
double V
Definition vpbe.h:134
double bulkIonicStrength
Definition vpbe.h:99
double L
Definition vpbe.h:132
int paramFlag
Definition vpbe.h:125
int ipkey
Definition vpbe.h:122
double ionRadii[MAXION]
Definition vpbe.h:105
double soluteCharge
Definition vpbe.h:118
double xkappa
Definition vpbe.h:108
double deblen
Definition vpbe.h:109
Vacc * acc
Definition vpbe.h:90
double z_mem
Definition vpbe.h:131
double solventDiel
Definition vpbe.h:94
int numIon
Definition vpbe.h:103
double membraneDiel
Definition vpbe.h:133
Vclist * clist
Definition vpbe.h:89
double zkappa2
Definition vpbe.h:110
double soluteZlen
Definition vpbe.h:117
double soluteCenter[3]
Definition vpbe.h:113
double soluteRadius
Definition vpbe.h:114
double soluteDiel
Definition vpbe.h:93
double soluteXlen
Definition vpbe.h:115
double T
Definition vpbe.h:92
int param2Flag
Definition vpbe.h:135
double maxIonRadius
Definition vpbe.h:100
double solventRadius
Definition vpbe.h:95
Vmem * vmem
Definition vpbe.h:86
Valist * alist
Definition vpbe.h:88
double ionConc[MAXION]
Definition vpbe.h:104
double smvolume
Definition vpbe.h:120
double zmagic
Definition vpbe.h:111
double ionQ[MAXION]
Definition vpbe.h:106
Contains declarations for class Vpbe.