APBS  1.4.1
routines.c
Go to the documentation of this file.
1 
54 #include "routines.h"
55 
56 VEMBED(rcsid="$Id$")
57 
58 VPUBLIC void startVio() { Vio_start(); }
59 
60 VPUBLIC Vparam* loadParameter(NOsh *nosh) {
61 
62  Vparam *param = VNULL;
63 
64  if (nosh->gotparm) {
65  param = Vparam_ctor();
66  switch (nosh->parmfmt) {
67  case NPF_FLAT:
68  Vnm_tprint( 1, "Reading parameter data from %s.\n",
69  nosh->parmpath);
70  if (Vparam_readFlatFile(param, "FILE", "ASC", VNULL,
71  nosh->parmpath) != 1) {
72  Vnm_tprint(2, "Error reading parameter file (%s)!\n", nosh->parmpath);
73  return VNULL;
74  }
75  break;
76  case NPF_XML:
77  Vnm_tprint( 1, "Reading parameter data from %s.\n",
78  nosh->parmpath);
79  if (Vparam_readXMLFile(param, "FILE", "ASC", VNULL,
80  nosh->parmpath) != 1) {
81  Vnm_tprint(2, "Error reading parameter file (%s)!\n", nosh->parmpath);
82  return VNULL;
83  }
84  break;
85  default:
86  Vnm_tprint(2, "Error! Undefined parameter file type (%d)!\n", nosh->parmfmt);
87  return VNULL;
88  } /* switch parmfmt */
89  }
90 
91  return param;
92 }
93 
94 
95 VPUBLIC int loadMolecules(NOsh *nosh, Vparam *param, Valist *alist[NOSH_MAXMOL]) {
96 
97  int i;
98  int use_params = 0;
99  Vrc_Codes rc;
100 
101  Vio *sock = VNULL;
102 
103  Vnm_tprint( 1, "Got paths for %d molecules\n", nosh->nmol);
104  if (nosh->nmol <= 0) {
105  Vnm_tprint(2, "You didn't specify any molecules (correctly)!\n");
106  Vnm_tprint(2, "Bailing out!\n");
107  return 0;
108  }
109 
110  if (nosh->gotparm) {
111  if (param == VNULL) {
112  Vnm_tprint(2, "Error! You don't have a valid parameter object!\n");
113  return 0;
114  }
115  use_params = 1;
116  }
117 
118  for (i=0; i<nosh->nmol; i++) {
119  if(alist[i] == VNULL){
120  alist[i] = Valist_ctor();
121  }else{
122  alist[i] = VNULL;
123  alist[i] = Valist_ctor();
124  }
125 
126  switch (nosh->molfmt[i]) {
127  case NMF_PQR:
128  /* Print out a warning to the user letting them know that we are overriding PQR
129  values for charge, radius and epsilon */
130  if (use_params) {
131  Vnm_print(2, "\nWARNING!! Radius/charge information from PQR file %s\n", nosh->molpath[i]);
132  Vnm_print(2, "will be replaced with data from parameter file (%s)!\n", nosh->parmpath);
133  }
134  Vnm_tprint( 1, "Reading PQR-format atom data from %s.\n",
135  nosh->molpath[i]);
136  sock = Vio_ctor("FILE", "ASC", VNULL, nosh->molpath[i], "r");
137  if (sock == VNULL) {
138  Vnm_print(2, "Problem opening virtual socket %s!\n",
139  nosh->molpath[i]);
140  return 0;
141  }
142  if (Vio_accept(sock, 0) < 0) {
143  Vnm_print(2, "Problem accepting virtual socket %s!\n",
144  nosh->molpath[i]);
145  return 0;
146  }
147  if(use_params){
148  rc = Valist_readPQR(alist[i], param, sock);
149  }else{
150  rc = Valist_readPQR(alist[i], VNULL, sock);
151  }
152  if(rc == 0) return 0;
153 
154  Vio_acceptFree(sock);
155  Vio_dtor(&sock);
156  break;
157  case NMF_PDB:
158  /* Load parameters */
159  if (!nosh->gotparm) {
160  Vnm_tprint(2, "NOsh: Error! Can't read PDB without specifying PARM file!\n");
161  return 0;
162  }
163  Vnm_tprint( 1, "Reading PDB-format atom data from %s.\n",
164  nosh->molpath[i]);
165  sock = Vio_ctor("FILE", "ASC", VNULL, nosh->molpath[i], "r");
166  if (sock == VNULL) {
167  Vnm_print(2, "Problem opening virtual socket %s!\n",
168  nosh->molpath[i]);
169  return 0;
170  }
171  if (Vio_accept(sock, 0) < 0) {
172  Vnm_print(2, "Problem accepting virtual socket %s!\n", nosh->molpath[i]);
173  return 0;
174  }
175  rc = Valist_readPDB(alist[i], param, sock);
176  /* If we are looking for an atom/residue that does not exist
177  * then abort and return 0 */
178  if(rc == 0)
179  return 0;
180 
181  Vio_acceptFree(sock);
182  Vio_dtor(&sock);
183  break;
184  case NMF_XML:
185  Vnm_tprint( 1, "Reading XML-format atom data from %s.\n",
186  nosh->molpath[i]);
187  sock = Vio_ctor("FILE", "ASC", VNULL, nosh->molpath[i], "r");
188  if (sock == VNULL) {
189  Vnm_print(2, "Problem opening virtual socket %s!\n",
190  nosh->molpath[i]);
191  return 0;
192  }
193  if (Vio_accept(sock, 0) < 0) {
194  Vnm_print(2, "Problem accepting virtual socket %s!\n",
195  nosh->molpath[i]);
196  return 0;
197  }
198  if(use_params){
199  rc = Valist_readXML(alist[i], param, sock);
200  }else{
201  rc = Valist_readXML(alist[i], VNULL, sock);
202  }
203  if(rc == 0)
204  return 0;
205 
206  Vio_acceptFree(sock);
207  Vio_dtor(&sock);
208  break;
209  default:
210  Vnm_tprint(2, "NOsh: Error! Undefined molecule file type \
211 (%d)!\n", nosh->molfmt[i]);
212  return 0;
213  } /* switch molfmt */
214 
215  if (rc != 1) {
216  Vnm_tprint( 2, "Error while reading molecule from %s\n",
217  nosh->molpath[i]);
218  return 0;
219  }
220 
221  Vnm_tprint( 1, " %d atoms\n", Valist_getNumberAtoms(alist[i]));
222  Vnm_tprint( 1, " Centered at (%4.3e, %4.3e, %4.3e)\n",
223  alist[i]->center[0], alist[i]->center[1],
224  alist[i]->center[2]);
225  Vnm_tprint( 1, " Net charge %3.2e e\n", alist[i]->charge);
226 
227  }
228 
229  return 1;
230 
231 }
232 
233 VPUBLIC void killMolecules(NOsh *nosh, Valist *alist[NOSH_MAXMOL]) {
234 
235  int i;
236 
237 #ifndef VAPBSQUIET
238  Vnm_tprint( 1, "Destroying %d molecules\n", nosh->nmol);
239 #endif
240 
241  for (i=0; i<nosh->nmol; i++)
242  Valist_dtor(&(alist[i]));
243 
244 }
245 
250 VPUBLIC int loadDielMaps(NOsh *nosh,
251  Vgrid *dielXMap[NOSH_MAXMOL],
252  Vgrid *dielYMap[NOSH_MAXMOL],
253  Vgrid *dielZMap[NOSH_MAXMOL]
254  ) {
255 
256  int i,
257  ii,
258  nx,
259  ny,
260  nz;
261  double sum,
262  hx,
263  hy,
264  hzed,
265  xmin,
266  ymin,
267  zmin;
268 
269  // Check to be sure we have dieletric map paths; if not, return.
270  if (nosh->ndiel > 0)
271  Vnm_tprint( 1, "Got paths for %d dielectric map sets\n",
272  nosh->ndiel);
273  else
274  return 1;
275 
276  // For each dielectric map path, read the data and calculate needed values.
277  for (i=0; i<nosh->ndiel; i++) {
278  Vnm_tprint( 1, "Reading x-shifted dielectric map data from \
279 %s:\n", nosh->dielXpath[i]);
280  dielXMap[i] = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
281 
282  // Determine the format and read data if the format is valid.
283  switch (nosh->dielfmt[i]) {
284  // OpenDX (Data Explorer) format
285  case VDF_DX:
286  if (Vgrid_readDX(dielXMap[i], "FILE", "ASC", VNULL,
287  nosh->dielXpath[i]) != 1) {
288  Vnm_tprint( 2, "Fatal error while reading from %s\n",
289  nosh->dielXpath[i]);
290  return 0;
291  }
292 
293  // Set grid sizes
294  nx = dielXMap[i]->nx;
295  ny = dielXMap[i]->ny;
296  nz = dielXMap[i]->nz;
297 
298  // Set spacings
299  hx = dielXMap[i]->hx;
300  hy = dielXMap[i]->hy;
301  hzed = dielXMap[i]->hzed;
302 
303  // Set minimum lower corner
304  xmin = dielXMap[i]->xmin;
305  ymin = dielXMap[i]->ymin;
306  zmin = dielXMap[i]->zmin;
307  Vnm_tprint(1, " %d x %d x %d grid\n", nx, ny, nz);
308  Vnm_tprint(1, " (%g, %g, %g) A spacings\n", hx, hy, hzed);
309  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
310  xmin, ymin, zmin);
311  sum = 0;
312  for (ii=0; ii<(nx*ny*nz); ii++)
313  sum += (dielXMap[i]->data[ii]);
314  sum = sum*hx*hy*hzed;
315  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
316  break;
317  // Binary file (GZip)
318  case VDF_GZ:
319  if (Vgrid_readGZ(dielXMap[i], nosh->dielXpath[i]) != 1) {
320  Vnm_tprint( 2, "Fatal error while reading from %s\n",
321  nosh->dielXpath[i]);
322  return 0;
323  }
324 
325  // Set grid sizes
326  nx = dielXMap[i]->nx;
327  ny = dielXMap[i]->ny;
328  nz = dielXMap[i]->nz;
329 
330  // Set spacings
331  hx = dielXMap[i]->hx;
332  hy = dielXMap[i]->hy;
333  hzed = dielXMap[i]->hzed;
334 
335  // Set minimum lower corner
336  xmin = dielXMap[i]->xmin;
337  ymin = dielXMap[i]->ymin;
338  zmin = dielXMap[i]->zmin;
339  Vnm_tprint(1, " %d x %d x %d grid\n", nx, ny, nz);
340  Vnm_tprint(1, " (%g, %g, %g) A spacings\n", hx, hy, hzed);
341  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
342  xmin, ymin, zmin);
343  sum = 0;
344  for (ii=0; ii<(nx*ny*nz); ii++)
345  sum += (dielXMap[i]->data[ii]);
346  sum = sum*hx*hy*hzed;
347  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
348  break;
349  // UHBD format
350  case VDF_UHBD:
351  Vnm_tprint( 2, "UHBD input not supported yet!\n");
352  return 0;
353  // AVS UCD format
354  case VDF_AVS:
355  Vnm_tprint( 2, "AVS input not supported yet!\n");
356  return 0;
357  // FEtk MC Simplex Format (MCSF)
358  case VDF_MCSF:
359  Vnm_tprint( 2, "MCSF input not supported yet!\n");
360  return 0;
361  default:
362  Vnm_tprint( 2, "Invalid data format (%d)!\n",
363  nosh->dielfmt[i]);
364  return 0;
365  }
366  Vnm_tprint( 1, "Reading y-shifted dielectric map data from \
367 %s:\n", nosh->dielYpath[i]);
368  dielYMap[i] = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
369 
370  // Determine the format and read data if the format is valid.
371  switch (nosh->dielfmt[i]) {
372  // OpenDX (Data Explorer) format
373  case VDF_DX:
374  if (Vgrid_readDX(dielYMap[i], "FILE", "ASC", VNULL,
375  nosh->dielYpath[i]) != 1) {
376  Vnm_tprint( 2, "Fatal error while reading from %s\n",
377  nosh->dielYpath[i]);
378  return 0;
379  }
380 
381  // Read grid
382  nx = dielYMap[i]->nx;
383  ny = dielYMap[i]->ny;
384  nz = dielYMap[i]->nz;
385 
386  // Read spacings
387  hx = dielYMap[i]->hx;
388  hy = dielYMap[i]->hy;
389  hzed = dielYMap[i]->hzed;
390 
391  // Read minimum lower corner
392  xmin = dielYMap[i]->xmin;
393  ymin = dielYMap[i]->ymin;
394  zmin = dielYMap[i]->zmin;
395  Vnm_tprint(1, " %d x %d x %d grid\n", nx, ny, nz);
396  Vnm_tprint(1, " (%g, %g, %g) A spacings\n", hx, hy, hzed);
397  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
398  xmin, ymin, zmin);
399  sum = 0;
400  for (ii=0; ii<(nx*ny*nz); ii++)
401  sum += (dielYMap[i]->data[ii]);
402  sum = sum*hx*hy*hzed;
403  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
404  break;
405  // Binary file (GZip) format
406  case VDF_GZ:
407  if (Vgrid_readGZ(dielYMap[i], nosh->dielYpath[i]) != 1) {
408  Vnm_tprint( 2, "Fatal error while reading from %s\n",
409  nosh->dielYpath[i]);
410  return 0;
411  }
412 
413  // Read grid
414  nx = dielYMap[i]->nx;
415  ny = dielYMap[i]->ny;
416  nz = dielYMap[i]->nz;
417 
418  // Read spacings
419  hx = dielYMap[i]->hx;
420  hy = dielYMap[i]->hy;
421  hzed = dielYMap[i]->hzed;
422 
423  // Read minimum lower corner
424  xmin = dielYMap[i]->xmin;
425  ymin = dielYMap[i]->ymin;
426  zmin = dielYMap[i]->zmin;
427  Vnm_tprint(1, " %d x %d x %d grid\n", nx, ny, nz);
428  Vnm_tprint(1, " (%g, %g, %g) A spacings\n", hx, hy, hzed);
429  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
430  xmin, ymin, zmin);
431  sum = 0;
432  for (ii=0; ii<(nx*ny*nz); ii++)
433  sum += (dielYMap[i]->data[ii]);
434  sum = sum*hx*hy*hzed;
435  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
436  break;
437  // UHBD format
438  case VDF_UHBD:
439  Vnm_tprint( 2, "UHBD input not supported yet!\n");
440  return 0;
441  // AVS UCD format
442  case VDF_AVS:
443  Vnm_tprint( 2, "AVS input not supported yet!\n");
444  return 0;
445  // FEtk MC Simplex Format (MCSF)
446  case VDF_MCSF:
447  Vnm_tprint( 2, "MCSF input not supported yet!\n");
448  return 0;
449  default:
450  Vnm_tprint( 2, "Invalid data format (%d)!\n",
451  nosh->dielfmt[i]);
452  return 0;
453  }
454 
455  Vnm_tprint( 1, "Reading z-shifted dielectric map data from \
456 %s:\n", nosh->dielZpath[i]);
457  dielZMap[i] = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
458 
459  // Determine the format and read data if the format is valid.
460  switch (nosh->dielfmt[i]) {
461  // OpenDX (Data Explorer) format
462  case VDF_DX:
463  if (Vgrid_readDX(dielZMap[i], "FILE", "ASC", VNULL,
464  nosh->dielZpath[i]) != 1) {
465  Vnm_tprint( 2, "Fatal error while reading from %s\n",
466  nosh->dielZpath[i]);
467  return 0;
468  }
469 
470  // Read grid
471  nx = dielZMap[i]->nx;
472  ny = dielZMap[i]->ny;
473  nz = dielZMap[i]->nz;
474 
475  // Read spacings
476  hx = dielZMap[i]->hx;
477  hy = dielZMap[i]->hy;
478  hzed = dielZMap[i]->hzed;
479 
480  // Read minimum lower corner
481  xmin = dielZMap[i]->xmin;
482  ymin = dielZMap[i]->ymin;
483  zmin = dielZMap[i]->zmin;
484  Vnm_tprint(1, " %d x %d x %d grid\n",
485  nx, ny, nz);
486  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
487  hx, hy, hzed);
488  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
489  xmin, ymin, zmin);
490  sum = 0;
491  for (ii=0; ii<(nx*ny*nz); ii++) sum += (dielZMap[i]->data[ii]);
492  sum = sum*hx*hy*hzed;
493  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
494  break;
495  // Binary file (GZip) format
496  case VDF_GZ:
497  if (Vgrid_readGZ(dielZMap[i], nosh->dielZpath[i]) != 1) {
498  Vnm_tprint( 2, "Fatal error while reading from %s\n",
499  nosh->dielZpath[i]);
500  return 0;
501  }
502 
503  // Read grid
504  nx = dielZMap[i]->nx;
505  ny = dielZMap[i]->ny;
506  nz = dielZMap[i]->nz;
507 
508  // Read spacings
509  hx = dielZMap[i]->hx;
510  hy = dielZMap[i]->hy;
511  hzed = dielZMap[i]->hzed;
512 
513  // Read minimum lower corner
514  xmin = dielZMap[i]->xmin;
515  ymin = dielZMap[i]->ymin;
516  zmin = dielZMap[i]->zmin;
517  Vnm_tprint(1, " %d x %d x %d grid\n",
518  nx, ny, nz);
519  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
520  hx, hy, hzed);
521  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
522  xmin, ymin, zmin);
523  sum = 0;
524  for (ii=0; ii<(nx*ny*nz); ii++) sum += (dielZMap[i]->data[ii]);
525  sum = sum*hx*hy*hzed;
526  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
527  break;
528  // UHBD format
529  case VDF_UHBD:
530  Vnm_tprint( 2, "UHBD input not supported yet!\n");
531  return 0;
532  // AVS UCD format
533  case VDF_AVS:
534  Vnm_tprint( 2, "AVS input not supported yet!\n");
535  return 0;
536  // FEtk MC Simplex Format (MCSF)
537  case VDF_MCSF:
538  Vnm_tprint( 2, "MCSF input not supported yet!\n");
539  return 0;
540  default:
541  Vnm_tprint( 2, "Invalid data format (%d)!\n",
542  nosh->dielfmt[i]);
543  return 0;
544  }
545  }
546 
547  return 1;
548 }
549 
550 VPUBLIC void killDielMaps(NOsh *nosh,
551  Vgrid *dielXMap[NOSH_MAXMOL],
552  Vgrid *dielYMap[NOSH_MAXMOL],
553  Vgrid *dielZMap[NOSH_MAXMOL]) {
554 
555  int i;
556 
557  if (nosh->ndiel > 0) {
558 #ifndef VAPBSQUIET
559  Vnm_tprint( 1, "Destroying %d dielectric map sets\n",
560  nosh->ndiel);
561 #endif
562  for (i=0; i<nosh->ndiel; i++) {
563  Vgrid_dtor(&(dielXMap[i]));
564  Vgrid_dtor(&(dielYMap[i]));
565  Vgrid_dtor(&(dielZMap[i]));
566  }
567  }
568  else return;
569 
570 }
571 
575 VPUBLIC int loadKappaMaps(NOsh *nosh,
576  Vgrid *map[NOSH_MAXMOL]
577  ) {
578 
579  int i,
580  ii,
581  len;
582  double sum;
583 
584  if (nosh->nkappa > 0)
585  Vnm_tprint( 1, "Got paths for %d kappa maps\n", nosh->nkappa);
586  else return 1;
587 
588  for (i=0; i<nosh->nkappa; i++) {
589  Vnm_tprint( 1, "Reading kappa map data from %s:\n",
590  nosh->kappapath[i]);
591  map[i] = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
592 
593  // Determine the format and read data if the format is valid.
594  switch (nosh->kappafmt[i]) {
595  // OpenDX (Data Explorer) format
596  case VDF_DX:
597  if (Vgrid_readDX(map[i], "FILE", "ASC", VNULL,
598  nosh->kappapath[i]) != 1) {
599  Vnm_tprint( 2, "Fatal error while reading from %s\n",
600  nosh->kappapath[i]);
601  return 0;
602  }
603  Vnm_tprint(1, " %d x %d x %d grid\n",
604  map[i]->nx, map[i]->ny, map[i]->nz);
605  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
606  map[i]->hx, map[i]->hy, map[i]->hzed);
607  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
608  map[i]->xmin, map[i]->ymin, map[i]->zmin);
609  sum = 0;
610  for (ii = 0, len = map[i]->nx * map[i]->ny * map[i]->nz;
611  ii < len;
612  ii++
613  ) {
614  sum += (map[i]->data[ii]);
615  }
616  sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
617  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
618  break;
619  // UHBD format
620  case VDF_UHBD:
621  Vnm_tprint( 2, "UHBD input not supported yet!\n");
622  return 0;
623  // FEtk MC Simplex Format (MCSF)
624  case VDF_MCSF:
625  Vnm_tprint( 2, "MCSF input not supported yet!\n");
626  return 0;
627  // AVS UCD format
628  case VDF_AVS:
629  Vnm_tprint( 2, "AVS input not supported yet!\n");
630  return 0;
631  // Binary file (GZip) format
632  case VDF_GZ:
633  if (Vgrid_readGZ(map[i], nosh->kappapath[i]) != 1) {
634  Vnm_tprint( 2, "Fatal error while reading from %s\n",
635  nosh->kappapath[i]);
636  return 0;
637  }
638  Vnm_tprint(1, " %d x %d x %d grid\n",
639  map[i]->nx, map[i]->ny, map[i]->nz);
640  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
641  map[i]->hx, map[i]->hy, map[i]->hzed);
642  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
643  map[i]->xmin, map[i]->ymin, map[i]->zmin);
644  sum = 0;
645  for (ii=0, len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
646  sum += (map[i]->data[ii]);
647  }
648  sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
649  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
650  break;
651  default:
652  Vnm_tprint( 2, "Invalid data format (%d)!\n",
653  nosh->kappafmt[i]);
654  return 0;
655  }
656  }
657 
658  return 1;
659 
660 }
661 
662 VPUBLIC void killKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL]) {
663 
664  int i;
665 
666  if (nosh->nkappa > 0) {
667 #ifndef VAPBSQUIET
668  Vnm_tprint( 1, "Destroying %d kappa maps\n", nosh->nkappa);
669 #endif
670  for (i=0; i<nosh->nkappa; i++) Vgrid_dtor(&(map[i]));
671  }
672  else return;
673 
674 }
675 
679 VPUBLIC int loadPotMaps(NOsh *nosh,
680  Vgrid *map[NOSH_MAXMOL]
681  ) {
682 
683  int i,
684  ii,
685  len;
686  double sum;
687 
688  if (nosh->npot > 0)
689  Vnm_tprint( 1, "Got paths for %d potential maps\n", nosh->npot);
690  else return 1;
691 
692  for (i=0; i<nosh->npot; i++) {
693  Vnm_tprint( 1, "Reading potential map data from %s:\n",
694  nosh->potpath[i]);
695  map[i] = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
696  switch (nosh->potfmt[i]) {
697  // OpenDX (Data Explorer) format
698  case VDF_DX:
699  // Binary file (GZip) format
700  case VDF_GZ:
701  if (nosh->potfmt[i] == VDF_DX) {
702  if (Vgrid_readDX(map[i], "FILE", "ASC", VNULL,
703  nosh->potpath[i]) != 1) {
704  Vnm_tprint( 2, "Fatal error while reading from %s\n",
705  nosh->potpath[i]);
706  return 0;
707  }
708  }else {
709  if (Vgrid_readGZ(map[i], nosh->potpath[i]) != 1) {
710  Vnm_tprint( 2, "Fatal error while reading from %s\n",
711  nosh->potpath[i]);
712  return 0;
713  }
714  }
715  Vnm_tprint(1, " %d x %d x %d grid\n",
716  map[i]->nx, map[i]->ny, map[i]->nz);
717  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
718  map[i]->hx, map[i]->hy, map[i]->hzed);
719  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
720  map[i]->xmin, map[i]->ymin, map[i]->zmin);
721  sum = 0;
722  for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
723  sum += (map[i]->data[ii]);
724  }
725  sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
726  Vnm_tprint(1, " Volume integral = %3.2e A^3\n", sum);
727  break;
728  // UHBD format
729  case VDF_UHBD:
730  Vnm_tprint( 2, "UHBD input not supported yet!\n");
731  return 0;
732  // FEtk MC Simplex Format (MCSF)
733  case VDF_MCSF:
734  Vnm_tprint( 2, "MCSF input not supported yet!\n");
735  return 0;
736  // AVS UCD format
737  case VDF_AVS:
738  Vnm_tprint( 2, "AVS input not supported yet!\n");
739  return 0;
740  default:
741  Vnm_tprint( 2, "Invalid data format (%d)!\n",
742  nosh->potfmt[i]);
743  return 0;
744  }
745  }
746 
747  return 1;
748 
749 }
750 
751 VPUBLIC void killPotMaps(NOsh *nosh,
752  Vgrid *map[NOSH_MAXMOL]
753  ) {
754 
755  int i;
756 
757  if (nosh->npot > 0) {
758 #ifndef VAPBSQUIET
759  Vnm_tprint( 1, "Destroying %d potential maps\n", nosh->npot);
760 #endif
761  for (i=0; i<nosh->npot; i++) Vgrid_dtor(&(map[i]));
762  }
763  else return;
764 
765 }
766 
770 VPUBLIC int loadChargeMaps(NOsh *nosh,
771  Vgrid *map[NOSH_MAXMOL]
772  ) {
773 
774  int i,
775  ii,
776  len;
777  double sum;
778 
779  if (nosh->ncharge > 0)
780  Vnm_tprint( 1, "Got paths for %d charge maps\n", nosh->ncharge);
781  else return 1;
782 
783  for (i=0; i<nosh->ncharge; i++) {
784  Vnm_tprint( 1, "Reading charge map data from %s:\n",
785  nosh->chargepath[i]);
786  map[i] = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, VNULL);
787 
788  // Determine data format and read data
789  switch (nosh->chargefmt[i]) {
790  case VDF_DX:
791  if (Vgrid_readDX(map[i], "FILE", "ASC", VNULL,
792  nosh->chargepath[i]) != 1) {
793  Vnm_tprint( 2, "Fatal error while reading from %s\n",
794  nosh->chargepath[i]);
795  return 0;
796  }
797  Vnm_tprint(1, " %d x %d x %d grid\n",
798  map[i]->nx, map[i]->ny, map[i]->nz);
799  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
800  map[i]->hx, map[i]->hy, map[i]->hzed);
801  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
802  map[i]->xmin, map[i]->ymin, map[i]->zmin);
803  sum = 0;
804  for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
805  sum += (map[i]->data[ii]);
806  }
807  sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
808  Vnm_tprint(1, " Charge map integral = %3.2e e\n", sum);
809  break;
810  case VDF_UHBD:
811  Vnm_tprint( 2, "UHBD input not supported yet!\n");
812  return 0;
813  case VDF_AVS:
814  Vnm_tprint( 2, "AVS input not supported yet!\n");
815  return 0;
816  case VDF_MCSF:
817  Vnm_tprint(2, "MCSF input not supported yet!\n");
818  return 0;
819  case VDF_GZ:
820  if (Vgrid_readGZ(map[i], nosh->chargepath[i]) != 1) {
821  Vnm_tprint( 2, "Fatal error while reading from %s\n",
822  nosh->chargepath[i]);
823  return 0;
824  }
825  Vnm_tprint(1, " %d x %d x %d grid\n",
826  map[i]->nx, map[i]->ny, map[i]->nz);
827  Vnm_tprint(1, " (%g, %g, %g) A spacings\n",
828  map[i]->hx, map[i]->hy, map[i]->hzed);
829  Vnm_tprint(1, " (%g, %g, %g) A lower corner\n",
830  map[i]->xmin, map[i]->ymin, map[i]->zmin);
831  sum = 0;
832  for (ii=0,len=map[i]->nx*map[i]->ny*map[i]->nz; ii<len; ii++) {
833  sum += (map[i]->data[ii]);
834  }
835  sum = sum*map[i]->hx*map[i]->hy*map[i]->hzed;
836  Vnm_tprint(1, " Charge map integral = %3.2e e\n", sum);
837  break;
838  default:
839  Vnm_tprint( 2, "Invalid data format (%d)!\n",
840  nosh->kappafmt[i]);
841  return 0;
842  }
843  }
844 
845  return 1;
846 
847 }
848 
849 VPUBLIC void killChargeMaps(NOsh *nosh,
850  Vgrid *map[NOSH_MAXMOL]
851  ) {
852 
853  int i;
854 
855  if (nosh->ncharge > 0) {
856 #ifndef VAPBSQUIET
857  Vnm_tprint( 1, "Destroying %d charge maps\n", nosh->ncharge);
858 #endif
859 
860  for (i=0; i<nosh->ncharge; i++) Vgrid_dtor(&(map[i]));
861  }
862 
863  else return;
864 
865 }
866 
867 VPUBLIC void printPBEPARM(PBEparm *pbeparm) {
868 
869  int i;
870  double ionstr = 0.0;
871 
872  for (i=0; i<pbeparm->nion; i++)
873  ionstr += 0.5*(VSQR(pbeparm->ionq[i])*pbeparm->ionc[i]);
874 
875  Vnm_tprint( 1, " Molecule ID: %d\n", pbeparm->molid);
876  switch (pbeparm->pbetype) {
877  case PBE_NPBE:
878  Vnm_tprint( 1, " Nonlinear traditional PBE\n");
879  break;
880  case PBE_LPBE:
881  Vnm_tprint( 1, " Linearized traditional PBE\n");
882  break;
883  case PBE_NRPBE:
884  Vnm_tprint( 1, " Nonlinear regularized PBE\n");
885  Vnm_tprint( 2, " ** Sorry, but Nathan broke the nonlinear regularized PBE implementation. **\n");
886  Vnm_tprint( 2, " ** Please let us know if you are interested in using it. **\n");
887  VASSERT(0);
888  break;
889  case PBE_LRPBE:
890  Vnm_tprint( 1, " Linearized regularized PBE\n");
891  break;
892  case PBE_SMPBE: /* SMPBE Added */
893  Vnm_tprint( 1, " Nonlinear Size-Modified PBE\n");
894  break;
895  default:
896  Vnm_tprint(2, " Unknown PBE type (%d)!\n", pbeparm->pbetype);
897  break;
898  }
899  if (pbeparm->bcfl == BCFL_ZERO) {
900  Vnm_tprint( 1, " Zero boundary conditions\n");
901  } else if (pbeparm->bcfl == BCFL_SDH) {
902  Vnm_tprint( 1, " Single Debye-Huckel sphere boundary \
903 conditions\n");
904  } else if (pbeparm->bcfl == BCFL_MDH) {
905  Vnm_tprint( 1, " Multiple Debye-Huckel sphere boundary \
906 conditions\n");
907  } else if (pbeparm->bcfl == BCFL_FOCUS) {
908  Vnm_tprint( 1, " Boundary conditions from focusing\n");
909  } else if (pbeparm->bcfl == BCFL_MAP) {
910  Vnm_tprint( 1, " Boundary conditions from potential map\n");
911  } else if (pbeparm->bcfl == BCFL_MEM) {
912  Vnm_tprint( 1, " Membrane potential boundary conditions.\n");
913  }
914  Vnm_tprint( 1, " %d ion species (%4.3f M ionic strength):\n",
915  pbeparm->nion, ionstr);
916  for (i=0; i<pbeparm->nion; i++) {
917  Vnm_tprint( 1, " %4.3f A-radius, %4.3f e-charge, \
918 %4.3f M concentration\n",
919  pbeparm->ionr[i], pbeparm->ionq[i], pbeparm->ionc[i]);
920  }
921 
922  if(pbeparm->pbetype == PBE_SMPBE){ /* SMPBE Added */
923  Vnm_tprint( 1, " Lattice spacing: %4.3f A (SMPBE) \n", pbeparm->smvolume);
924  Vnm_tprint( 1, " Relative size parameter: %4.3f (SMPBE) \n", pbeparm->smsize);
925  }
926 
927  Vnm_tprint( 1, " Solute dielectric: %4.3f\n", pbeparm->pdie);
928  Vnm_tprint( 1, " Solvent dielectric: %4.3f\n", pbeparm->sdie);
929  switch (pbeparm->srfm) {
930  case 0:
931  Vnm_tprint( 1, " Using \"molecular\" surface \
932 definition; no smoothing\n");
933  Vnm_tprint( 1, " Solvent probe radius: %4.3f A\n",
934  pbeparm->srad);
935  break;
936  case 1:
937  Vnm_tprint( 1, " Using \"molecular\" surface definition;\
938 harmonic average smoothing\n");
939  Vnm_tprint( 1, " Solvent probe radius: %4.3f A\n",
940  pbeparm->srad);
941  break;
942  case 2:
943  Vnm_tprint( 1, " Using spline-based surface definition;\
944 window = %4.3f\n", pbeparm->swin);
945  break;
946  default:
947  break;
948  }
949  Vnm_tprint( 1, " Temperature: %4.3f K\n", pbeparm->temp);
950  if (pbeparm->calcenergy != PCE_NO) Vnm_tprint( 1, " Electrostatic \
951 energies will be calculated\n");
952  if (pbeparm->calcforce == PCF_TOTAL) Vnm_tprint( 1, " Net solvent \
953 forces will be calculated \n");
954  if (pbeparm->calcforce == PCF_COMPS) Vnm_tprint( 1, " All-atom \
955 solvent forces will be calculated\n");
956  for (i=0; i<pbeparm->numwrite; i++) {
957  switch (pbeparm->writetype[i]) {
958  case VDT_CHARGE:
959  Vnm_tprint(1, " Charge distribution to be written to ");
960  break;
961  case VDT_POT:
962  Vnm_tprint(1, " Potential to be written to ");
963  break;
964  case VDT_SMOL:
965  Vnm_tprint(1, " Molecular solvent accessibility \
966 to be written to ");
967  break;
968  case VDT_SSPL:
969  Vnm_tprint(1, " Spline-based solvent accessibility \
970 to be written to ");
971  break;
972  case VDT_VDW:
973  Vnm_tprint(1, " van der Waals solvent accessibility \
974 to be written to ");
975  break;
976  case VDT_IVDW:
977  Vnm_tprint(1, " Ion accessibility to be written to ");
978  break;
979  case VDT_LAP:
980  Vnm_tprint(1, " Potential Laplacian to be written to ");
981  break;
982  case VDT_EDENS:
983  Vnm_tprint(1, " Energy density to be written to ");
984  break;
985  case VDT_NDENS:
986  Vnm_tprint(1, " Ion number density to be written to ");
987  break;
988  case VDT_QDENS:
989  Vnm_tprint(1, " Ion charge density to be written to ");
990  break;
991  case VDT_DIELX:
992  Vnm_tprint(1, " X-shifted dielectric map to be written \
993 to ");
994  break;
995  case VDT_DIELY:
996  Vnm_tprint(1, " Y-shifted dielectric map to be written \
997 to ");
998  break;
999  case VDT_DIELZ:
1000  Vnm_tprint(1, " Z-shifted dielectric map to be written \
1001 to ");
1002  break;
1003  case VDT_KAPPA:
1004  Vnm_tprint(1, " Kappa map to be written to ");
1005  break;
1006  case VDT_ATOMPOT:
1007  Vnm_tprint(1, " Atom potentials to be written to ");
1008  break;
1009  default:
1010  Vnm_tprint(2, " Invalid data type for writing!\n");
1011  break;
1012  }
1013  switch (pbeparm->writefmt[i]) {
1014  case VDF_DX:
1015  Vnm_tprint(1, "%s.%s\n", pbeparm->writestem[i], "dx");
1016  break;
1017  case VDF_GZ:
1018  Vnm_tprint(1, "%s.%s\n", pbeparm->writestem[i], "dx.gz");
1019  break;
1020  case VDF_UHBD:
1021  Vnm_tprint(1, "%s.%s\n", pbeparm->writestem[i], "grd");
1022  break;
1023  case VDF_AVS:
1024  Vnm_tprint(1, "%s.%s\n", pbeparm->writestem[i], "ucd");
1025  break;
1026  case VDF_MCSF:
1027  Vnm_tprint(1, "%s.%s\n", pbeparm->writestem[i], "mcsf");
1028  break;
1029  case VDF_FLAT:
1030  Vnm_tprint(1, "%s.%s\n", pbeparm->writestem[i], "txt");
1031  break;
1032  default:
1033  Vnm_tprint(2, " Invalid format for writing!\n");
1034  break;
1035  }
1036 
1037  }
1038 
1039 }
1040 
1041 VPUBLIC void printMGPARM(MGparm *mgparm, double realCenter[3]) {
1042 
1043  switch (mgparm->chgm) {
1044  case 0:
1045  Vnm_tprint(1, " Using linear spline charge discretization.\n");
1046  break;
1047  case 1:
1048  Vnm_tprint(1, " Using cubic spline charge discretization.\n");
1049  break;
1050  default:
1051  break;
1052  }
1053  if (mgparm->type == MCT_PARALLEL) {
1054  Vnm_tprint( 1, " Partition overlap fraction = %g\n",
1055  mgparm->ofrac);
1056  Vnm_tprint( 1, " Processor array = %d x %d x %d\n",
1057  mgparm->pdime[0], mgparm->pdime[1], mgparm->pdime[2]);
1058  }
1059  Vnm_tprint( 1, " Grid dimensions: %d x %d x %d\n",
1060  mgparm->dime[0], mgparm->dime[1], mgparm->dime[2]);
1061  Vnm_tprint( 1, " Grid spacings: %4.3f x %4.3f x %4.3f\n",
1062  mgparm->grid[0], mgparm->grid[1], mgparm->grid[2]);
1063  Vnm_tprint( 1, " Grid lengths: %4.3f x %4.3f x %4.3f\n",
1064  mgparm->glen[0], mgparm->glen[1], mgparm->glen[2]);
1065  Vnm_tprint( 1, " Grid center: (%4.3f, %4.3f, %4.3f)\n",
1066  realCenter[0], realCenter[1], realCenter[2]);
1067  Vnm_tprint( 1, " Multigrid levels: %d\n", mgparm->nlev);
1068 
1069 }
1070 
1074 VPUBLIC int initMG(int icalc,
1075  NOsh *nosh, MGparm *mgparm,
1076  PBEparm *pbeparm,
1077  double realCenter[3],
1078  Vpbe *pbe[NOSH_MAXCALC],
1079  Valist *alist[NOSH_MAXMOL],
1080  Vgrid *dielXMap[NOSH_MAXMOL],
1081  Vgrid *dielYMap[NOSH_MAXMOL],
1082  Vgrid *dielZMap[NOSH_MAXMOL],
1083  Vgrid *kappaMap[NOSH_MAXMOL],
1084  Vgrid *chargeMap[NOSH_MAXMOL],
1085  Vpmgp *pmgp[NOSH_MAXCALC],
1086  Vpmg *pmg[NOSH_MAXCALC],
1087  Vgrid *potMap[NOSH_MAXMOL]
1088  ) {
1089 
1090  int j,
1091  focusFlag,
1092  iatom;
1093  size_t bytesTotal,
1094  highWater;
1095  double sparm,
1096  iparm,
1097  q;
1098  Vatom *atom = VNULL;
1099  Vgrid *theDielXMap = VNULL,
1100  *theDielYMap = VNULL,
1101  *theDielZMap = VNULL;
1102  Vgrid *theKappaMap = VNULL,
1103  *thePotMap = VNULL,
1104  *theChargeMap = VNULL;
1105  Valist *myalist = VNULL;
1106 
1107  Vnm_tstart(APBS_TIMER_SETUP, "Setup timer");
1108 
1109  /* Update the grid center */
1110  for (j=0; j<3; j++) realCenter[j] = mgparm->center[j];
1111 
1112  /* Check for completely-neutral molecule */
1113  q = 0;
1114  myalist = alist[pbeparm->molid-1];
1115  for (iatom=0; iatom<Valist_getNumberAtoms(myalist); iatom++) {
1116  atom = Valist_getAtom(myalist, iatom);
1117  q += VSQR(Vatom_getCharge(atom));
1118  }
1119  /* D. Gohara 10/22/09 - disabled
1120  if (q < (1e-6)) {
1121  Vnm_tprint(2, "Molecule #%d is uncharged!\n", pbeparm->molid);
1122  Vnm_tprint(2, "Sum square charge = %g!\n", q);
1123  return 0;
1124  }
1125  */
1126 
1127  /* Set up PBE object */
1128  Vnm_tprint(0, "Setting up PBE object...\n");
1129  if (pbeparm->srfm == VSM_SPLINE) {
1130  sparm = pbeparm->swin;
1131  } else {
1132  sparm = pbeparm->srad;
1133  }
1134  if (pbeparm->nion > 0) {
1135  iparm = pbeparm->ionr[0];
1136  } else {
1137  iparm = 0.0;
1138  }
1139  if (pbeparm->bcfl == BCFL_FOCUS) {
1140  if (icalc == 0) {
1141  Vnm_tprint( 2, "Can't focus first calculation!\n");
1142  return 0;
1143  }
1144  focusFlag = 1;
1145  } else {
1146  focusFlag = 0;
1147  }
1148 
1149  // Construct Vpbe object
1150  pbe[icalc] = Vpbe_ctor(myalist, pbeparm->nion,
1151  pbeparm->ionc, pbeparm->ionr, pbeparm->ionq,
1152  pbeparm->temp, pbeparm->pdie,
1153  pbeparm->sdie, sparm, focusFlag, pbeparm->sdens,
1154  pbeparm->zmem, pbeparm->Lmem, pbeparm->mdie,
1155  pbeparm->memv);
1156 
1157  /* Set up PDE object */
1158  Vnm_tprint(0, "Setting up PDE object...\n");
1159  switch (pbeparm->pbetype) {
1160  case PBE_NPBE:
1161  /* TEMPORARY USEAQUA */
1162  mgparm->nonlintype = NONLIN_NPBE;
1163  mgparm->method = (mgparm->useAqua == 1) ? VSOL_NewtonAqua : VSOL_Newton;
1164  pmgp[icalc] = Vpmgp_ctor(mgparm);
1165  break;
1166  case PBE_LPBE:
1167  /* TEMPORARY USEAQUA */
1168  mgparm->nonlintype = NONLIN_LPBE;
1169  mgparm->method = (mgparm->useAqua == 1) ? VSOL_CGMGAqua : VSOL_MG;
1170  pmgp[icalc] = Vpmgp_ctor(mgparm);
1171  break;
1172  case PBE_LRPBE:
1173  Vnm_tprint(2, "Sorry, LRPBE isn't supported with the MG solver!\n");
1174  return 0;
1175  case PBE_NRPBE:
1176  Vnm_tprint(2, "Sorry, NRPBE isn't supported with the MG solver!\n");
1177  return 0;
1178  case PBE_SMPBE: /* SMPBE Added */
1179  mgparm->nonlintype = NONLIN_SMPBE;
1180  pmgp[icalc] = Vpmgp_ctor(mgparm);
1181 
1182  /* Copy Code */
1183  pbe[icalc]->smsize = pbeparm->smsize;
1184  pbe[icalc]->smvolume = pbeparm->smvolume;
1185  pbe[icalc]->ipkey = pmgp[icalc]->ipkey;
1186 
1187  break;
1188  default:
1189  Vnm_tprint(2, "Error! Unknown PBE type (%d)!\n", pbeparm->pbetype);
1190  return 0;
1191  }
1192  Vnm_tprint(0, "Setting PDE center to local center...\n");
1193  pmgp[icalc]->bcfl = pbeparm->bcfl;
1194  pmgp[icalc]->xcent = realCenter[0];
1195  pmgp[icalc]->ycent = realCenter[1];
1196  pmgp[icalc]->zcent = realCenter[2];
1197 
1198  if (pbeparm->bcfl == BCFL_FOCUS) {
1199  if (icalc == 0) {
1200  Vnm_tprint( 2, "Can't focus first calculation!\n");
1201  return 0;
1202  }
1203  /* Focusing requires the previous calculation in order to setup the
1204  current run... */
1205  pmg[icalc] = Vpmg_ctor(pmgp[icalc], pbe[icalc], 1, pmg[icalc-1],
1206  mgparm, pbeparm->calcenergy);
1207  /* ...however, it should be done with the previous calculation now, so
1208  we should be able to destroy it here. */
1209  /* Vpmg_dtor(&(pmg[icalc-1])); */
1210  } else {
1211  if (icalc>0) Vpmg_dtor(&(pmg[icalc-1]));
1212  pmg[icalc] = Vpmg_ctor(pmgp[icalc], pbe[icalc], 0, VNULL, mgparm, PCE_NO);
1213  }
1214  if (icalc>0) {
1215  Vpmgp_dtor(&(pmgp[icalc-1]));
1216  Vpbe_dtor(&(pbe[icalc-1]));
1217  }
1218  if (pbeparm->useDielMap) {
1219  if ((pbeparm->dielMapID-1) < nosh->ndiel) {
1220  theDielXMap = dielXMap[pbeparm->dielMapID-1];
1221  } else {
1222  Vnm_print(2, "Error! %d is not a valid dielectric map ID!\n",
1223  pbeparm->dielMapID);
1224  return 0;
1225  }
1226  }
1227  if (pbeparm->useDielMap) {
1228  if ((pbeparm->dielMapID-1) < nosh->ndiel) {
1229  theDielYMap = dielYMap[pbeparm->dielMapID-1];
1230  } else {
1231  Vnm_print(2, "Error! %d is not a valid dielectric map ID!\n",
1232  pbeparm->dielMapID);
1233  return 0;
1234  }
1235  }
1236  if (pbeparm->useDielMap) {
1237  if ((pbeparm->dielMapID-1) < nosh->ndiel) {
1238  theDielZMap = dielZMap[pbeparm->dielMapID-1];
1239  } else {
1240  Vnm_print(2, "Error! %d is not a valid dielectric map ID!\n",
1241  pbeparm->dielMapID);
1242  return 0;
1243  }
1244  }
1245  if (pbeparm->useKappaMap) {
1246  if ((pbeparm->kappaMapID-1) < nosh->nkappa) {
1247  theKappaMap = kappaMap[pbeparm->kappaMapID-1];
1248  } else {
1249  Vnm_print(2, "Error! %d is not a valid kappa map ID!\n",
1250  pbeparm->kappaMapID);
1251  return 0;
1252  }
1253  }
1254  if (pbeparm->usePotMap) {
1255  if ((pbeparm->potMapID-1) < nosh->npot) {
1256  thePotMap = potMap[pbeparm->potMapID-1];
1257  } else {
1258  Vnm_print(2, "Error! %d is not a valid potential map ID!\n",
1259  pbeparm->potMapID);
1260  return 0;
1261  }
1262  }
1263  if (pbeparm->useChargeMap) {
1264  if ((pbeparm->chargeMapID-1) < nosh->ncharge) {
1265  theChargeMap = chargeMap[pbeparm->chargeMapID-1];
1266  } else {
1267  Vnm_print(2, "Error! %d is not a valid charge map ID!\n",
1268  pbeparm->chargeMapID);
1269  return 0;
1270  }
1271  }
1272 
1273  if (pbeparm->bcfl == BCFL_MAP && thePotMap == VNULL) {
1274  Vnm_print(2, "Warning: You specified 'bcfl map' in the input file, but no potential map was found.\n");
1275  Vnm_print(2, " You must specify 'usemap pot' statement in the APBS input file!\n");
1276  Vnm_print(2, "Bailing out ...\n");
1277  return 0;
1278  }
1279 
1280  // Initialize calculation coefficients
1281  if (!Vpmg_fillco(pmg[icalc],
1282  pbeparm->srfm, pbeparm->swin, mgparm->chgm,
1283  pbeparm->useDielMap, theDielXMap,
1284  pbeparm->useDielMap, theDielYMap,
1285  pbeparm->useDielMap, theDielZMap,
1286  pbeparm->useKappaMap, theKappaMap,
1287  pbeparm->usePotMap, thePotMap,
1288  pbeparm->useChargeMap, theChargeMap)) {
1289  Vnm_print(2, "initMG: problems setting up coefficients (fillco)!\n");
1290  return 0;
1291  }
1292 
1293  /* Print a few derived parameters */
1294 #ifndef VAPBSQUIET
1295  Vnm_tprint(1, " Debye length: %g A\n", Vpbe_getDeblen(pbe[icalc]));
1296 #endif
1297 
1298  /* Setup time statistics */
1299  Vnm_tstop(APBS_TIMER_SETUP, "Setup timer");
1300 
1301  /* Memory statistics */
1302  bytesTotal = Vmem_bytesTotal();
1303  highWater = Vmem_highWaterTotal();
1304 
1305 #ifndef VAPBSQUIET
1306  Vnm_tprint( 1, " Current memory usage: %4.3f MB total, \
1307 %4.3f MB high water\n", (double)(bytesTotal)/(1024.*1024.),
1308  (double)(highWater)/(1024.*1024.));
1309 #endif
1310 
1311  return 1;
1312 
1313 }
1314 
1315 VPUBLIC void killMG(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC],
1316  Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC]) {
1317 
1318  int i;
1319 
1320 #ifndef VAPBSQUIET
1321  Vnm_tprint(1, "Destroying multigrid structures.\n");
1322 #endif
1323 
1324  /*
1325  There appears to be a relationship (or this is a bug in Linux, can't tell
1326  at the moment, since Linux is the only OS that seems to be affected)
1327  between one of the three object types: Vpbe, Vpmg or Vpmgp that requires
1328  deallocations to be performed in a specific order. This results in a
1329  bug some of the time when freeing Vpmg objects below. Therefore it
1330  appears to be important to release the Vpmg structs BEFORE the Vpmgp structs .
1331  */
1332  Vpmg_dtor(&(pmg[nosh->ncalc-1]));
1333 
1334  for(i=0;i<nosh->ncalc;i++){
1335  Vpbe_dtor(&(pbe[i]));
1336  Vpmgp_dtor(&(pmgp[i]));
1337  }
1338 
1339 }
1340 
1341 VPUBLIC int solveMG(NOsh *nosh,
1342  Vpmg *pmg,
1343  MGparm_CalcType type
1344  ) {
1345 
1346  int nx,
1347  ny,
1348  nz,
1349  i;
1350 
1351  if (nosh != VNULL) {
1352  if (nosh->bogus) return 1;
1353  }
1354 
1355  Vnm_tstart(APBS_TIMER_SOLVER, "Solver timer");
1356 
1357 
1358  if (type != MCT_DUMMY) {
1359  if (!Vpmg_solve(pmg)) {
1360  Vnm_print(2, " Error during PDE solution!\n");
1361  return 0;
1362  }
1363  } else {
1364  Vnm_tprint( 1," Skipping solve for mg-dummy run; zeroing \
1365 solution array\n");
1366  nx = pmg->pmgp->nx;
1367  ny = pmg->pmgp->ny;
1368  nz = pmg->pmgp->nz;
1369  for (i=0; i<nx*ny*nz; i++) pmg->u[i] = 0.0;
1370  }
1371  Vnm_tstop(APBS_TIMER_SOLVER, "Solver timer");
1372 
1373  return 1;
1374 
1375 }
1376 
1377 VPUBLIC int setPartMG(NOsh *nosh,
1378  MGparm *mgparm,
1379  Vpmg *pmg
1380  ) {
1381 
1382  int j;
1383  double partMin[3],
1384  partMax[3];
1385 
1386  if (nosh->bogus) return 1;
1387 
1388  if (mgparm->type == MCT_PARALLEL) {
1389  for (j=0; j<3; j++) {
1390  partMin[j] = mgparm->partDisjCenter[j] - 0.5*mgparm->partDisjLength[j];
1391  partMax[j] = mgparm->partDisjCenter[j] + 0.5*mgparm->partDisjLength[j];
1392  }
1393 #if 0
1394  Vnm_tprint(1, "setPartMG (%s, %d): Disj part center = (%g, %g, %g)\n",
1395  __FILE__, __LINE__,
1396  mgparm->partDisjCenter[0],
1397  mgparm->partDisjCenter[1],
1398  mgparm->partDisjCenter[2]
1399  );
1400  Vnm_tprint(1, "setPartMG (%s, %d): Disj part lower corner = (%g, %g, %g)\n",
1401  __FILE__, __LINE__, partMin[0], partMin[1], partMin[2]);
1402  Vnm_tprint(1, "setPartMG (%s, %d): Disj part upper corner = (%g, %g, %g)\n",
1403  __FILE__, __LINE__,
1404  partMax[0], partMax[1], partMax[2]);
1405 #endif
1406  } else {
1407  for (j=0; j<3; j++) {
1408  partMin[j] = mgparm->center[j] - 0.5*mgparm->glen[j];
1409  partMax[j] = mgparm->center[j] + 0.5*mgparm->glen[j];
1410  }
1411  }
1412  /* Vnm_print(1, "DEBUG (%s, %d): setPartMG calling setPart with upper corner \
1413 %g %g %g and lower corner %g %g %g\n", __FILE__, __LINE__,
1414  partMin[0], partMin[1], partMin[2],
1415  partMax[0], partMax[1], partMax[2]); */
1416  Vpmg_setPart(pmg, partMin, partMax, mgparm->partDisjOwnSide);
1417 
1418 
1419  return 1;
1420 
1421 }
1422 
1423 VPUBLIC int energyMG(NOsh *nosh,
1424  int icalc,
1425  Vpmg *pmg,
1426  int *nenergy,
1427  double *totEnergy,
1428  double *qfEnergy,
1429  double *qmEnergy,
1430  double *dielEnergy
1431  ) {
1432 
1433  Valist *alist;
1434  Vatom *atom;
1435  int i,
1436  extEnergy;
1437  double tenergy;
1438  MGparm *mgparm;
1439  PBEparm *pbeparm;
1440 
1441  mgparm = nosh->calc[icalc]->mgparm;
1442  pbeparm = nosh->calc[icalc]->pbeparm;
1443 
1444  Vnm_tstart(APBS_TIMER_ENERGY, "Energy timer");
1445  extEnergy = 1;
1446 
1447  if (pbeparm->calcenergy == PCE_TOTAL) {
1448  *nenergy = 1;
1449  /* Some processors don't count */
1450  if (nosh->bogus == 0) {
1451  *totEnergy = Vpmg_energy(pmg, extEnergy);
1452 #ifndef VAPBSQUIET
1453  Vnm_tprint( 1, " Total electrostatic energy = %1.12E kJ/mol\n",
1454  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*totEnergy));
1455 #endif
1456  } else *totEnergy = 0;
1457  } else if (pbeparm->calcenergy == PCE_COMPS) {
1458  *nenergy = 1;
1459  *totEnergy = Vpmg_energy(pmg, extEnergy);
1460  *qfEnergy = Vpmg_qfEnergy(pmg, extEnergy);
1461  *qmEnergy = Vpmg_qmEnergy(pmg, extEnergy);
1462  *dielEnergy = Vpmg_dielEnergy(pmg, extEnergy);
1463 #ifndef VAPBSQUIET
1464  Vnm_tprint( 1, " Total electrostatic energy = %1.12E \
1465 kJ/mol\n", Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*totEnergy));
1466  Vnm_tprint( 1, " Fixed charge energy = %g kJ/mol\n",
1467  0.5*Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*qfEnergy));
1468  Vnm_tprint( 1, " Mobile charge energy = %g kJ/mol\n",
1469  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*qmEnergy));
1470  Vnm_tprint( 1, " Dielectric energy = %g kJ/mol\n",
1471  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*dielEnergy));
1472  Vnm_tprint( 1, " Per-atom energies:\n");
1473 #endif
1474  alist = pmg->pbe->alist;
1475  for (i=0; i<Valist_getNumberAtoms(alist); i++) {
1476  atom = Valist_getAtom(alist, i);
1477  tenergy = Vpmg_qfAtomEnergy(pmg, atom);
1478 #ifndef VAPBSQUIET
1479  Vnm_tprint( 1, " Atom %d: %1.12E kJ/mol\n", i,
1480  0.5*Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*tenergy);
1481 #endif
1482  }
1483  } else *nenergy = 0;
1484 
1485  Vnm_tstop(APBS_TIMER_ENERGY, "Energy timer");
1486 
1487  return 1;
1488 }
1489 
1490 VPUBLIC int forceMG(Vmem *mem,
1491  NOsh *nosh,
1492  PBEparm *pbeparm,
1493  MGparm *mgparm,
1494  Vpmg *pmg,
1495  int *nforce,
1496  AtomForce **atomForce,
1497  Valist *alist[NOSH_MAXMOL]
1498  ) {
1499 
1500  int j,
1501  k;
1502  double qfForce[3],
1503  dbForce[3],
1504  ibForce[3];
1505 
1506  Vnm_tstart(APBS_TIMER_FORCE, "Force timer");
1507 
1508 #ifndef VAPBSQUIET
1509  Vnm_tprint( 1," Calculating forces...\n");
1510 #endif
1511 
1512  if (pbeparm->calcforce == PCF_TOTAL) {
1513  *nforce = 1;
1514  *atomForce = (AtomForce *)Vmem_malloc(mem, 1, sizeof(AtomForce));
1515  /* Clear out force arrays */
1516  for (j=0; j<3; j++) {
1517  (*atomForce)[0].qfForce[j] = 0;
1518  (*atomForce)[0].ibForce[j] = 0;
1519  (*atomForce)[0].dbForce[j] = 0;
1520  }
1521  for (j=0;j<Valist_getNumberAtoms(alist[pbeparm->molid-1]);j++) {
1522  if (nosh->bogus == 0) {
1523  VASSERT(Vpmg_qfForce(pmg, qfForce, j, mgparm->chgm));
1524  VASSERT(Vpmg_ibForce(pmg, ibForce, j, pbeparm->srfm));
1525  VASSERT(Vpmg_dbForce(pmg, dbForce, j, pbeparm->srfm));
1526  } else {
1527  for (k=0; k<3; k++) {
1528  qfForce[k] = 0;
1529  ibForce[k] = 0;
1530  dbForce[k] = 0;
1531  }
1532  }
1533  for (k=0; k<3; k++) {
1534  (*atomForce)[0].qfForce[k] += qfForce[k];
1535  (*atomForce)[0].ibForce[k] += ibForce[k];
1536  (*atomForce)[0].dbForce[k] += dbForce[k];
1537  }
1538  }
1539 #ifndef VAPBSQUIET
1540  Vnm_tprint( 1, " Printing net forces for molecule %d (kJ/mol/A)\n",
1541  pbeparm->molid);
1542  Vnm_tprint( 1, " Legend:\n");
1543  Vnm_tprint( 1, " qf -- fixed charge force\n");
1544  Vnm_tprint( 1, " db -- dielectric boundary force\n");
1545  Vnm_tprint( 1, " ib -- ionic boundary force\n");
1546  Vnm_tprint( 1, " qf %4.3e %4.3e %4.3e\n",
1547  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].qfForce[0],
1548  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].qfForce[1],
1549  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].qfForce[2]);
1550  Vnm_tprint( 1, " ib %4.3e %4.3e %4.3e\n",
1551  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].ibForce[0],
1552  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].ibForce[1],
1553  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].ibForce[2]);
1554  Vnm_tprint( 1, " db %4.3e %4.3e %4.3e\n",
1555  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].dbForce[0],
1556  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].dbForce[1],
1557  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*atomForce)[0].dbForce[2]);
1558 #endif
1559  } else if (pbeparm->calcforce == PCF_COMPS) {
1560  *nforce = Valist_getNumberAtoms(alist[pbeparm->molid-1]);
1561  *atomForce = (AtomForce *)Vmem_malloc(mem, *nforce,
1562  sizeof(AtomForce));
1563 #ifndef VAPBSQUIET
1564  Vnm_tprint( 1, " Printing per-atom forces for molecule %d (kJ/mol/A)\n",
1565  pbeparm->molid);
1566  Vnm_tprint( 1, " Legend:\n");
1567  Vnm_tprint( 1, " tot n -- total force for atom n\n");
1568  Vnm_tprint( 1, " qf n -- fixed charge force for atom n\n");
1569  Vnm_tprint( 1, " db n -- dielectric boundary force for atom n\n");
1570  Vnm_tprint( 1, " ib n -- ionic boundary force for atom n\n");
1571 #endif
1572  for (j=0;j<Valist_getNumberAtoms(alist[pbeparm->molid-1]);j++) {
1573  if (nosh->bogus == 0) {
1574  VASSERT(Vpmg_qfForce(pmg, (*atomForce)[j].qfForce, j,
1575  mgparm->chgm));
1576  VASSERT(Vpmg_ibForce(pmg, (*atomForce)[j].ibForce, j,
1577  pbeparm->srfm));
1578  VASSERT(Vpmg_dbForce(pmg, (*atomForce)[j].dbForce, j,
1579  pbeparm->srfm));
1580  } else {
1581  for (k=0; k<3; k++) {
1582  (*atomForce)[j].qfForce[k] = 0;
1583  (*atomForce)[j].ibForce[k] = 0;
1584  (*atomForce)[j].dbForce[k] = 0;
1585  }
1586  }
1587 #ifndef VAPBSQUIET
1588  Vnm_tprint( 1, "mgF tot %d %4.3e %4.3e %4.3e\n", j,
1589  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1590  *((*atomForce)[j].qfForce[0]+(*atomForce)[j].ibForce[0]+
1591  (*atomForce)[j].dbForce[0]),
1592  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1593  *((*atomForce)[j].qfForce[1]+(*atomForce)[j].ibForce[1]+
1594  (*atomForce)[j].dbForce[1]),
1595  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1596  *((*atomForce)[j].qfForce[2]+(*atomForce)[j].ibForce[2]+
1597  (*atomForce)[j].dbForce[2]));
1598  Vnm_tprint( 1, "mgF qf %d %4.3e %4.3e %4.3e\n", j,
1599  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1600  *(*atomForce)[j].qfForce[0],
1601  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1602  *(*atomForce)[j].qfForce[1],
1603  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1604  *(*atomForce)[j].qfForce[2]);
1605  Vnm_tprint( 1, "mgF ib %d %4.3e %4.3e %4.3e\n", j,
1606  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1607  *(*atomForce)[j].ibForce[0],
1608  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1609  *(*atomForce)[j].ibForce[1],
1610  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1611  *(*atomForce)[j].ibForce[2]);
1612  Vnm_tprint( 1, "mgF db %d %4.3e %4.3e %4.3e\n", j,
1613  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1614  *(*atomForce)[j].dbForce[0],
1615  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1616  *(*atomForce)[j].dbForce[1],
1617  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na \
1618  *(*atomForce)[j].dbForce[2]);
1619 #endif
1620  }
1621  } else *nforce = 0;
1622 
1623  Vnm_tstop(APBS_TIMER_FORCE, "Force timer");
1624 
1625  return 1;
1626 }
1627 
1628 VPUBLIC void killEnergy() {
1629 
1630 #ifndef VAPBSQUIET
1631  Vnm_tprint(1, "No energy arrays to destroy.\n");
1632 #endif
1633 
1634 }
1635 
1636 VPUBLIC void killForce(Vmem *mem, NOsh *nosh, int nforce[NOSH_MAXCALC],
1637  AtomForce *atomForce[NOSH_MAXCALC]) {
1638 
1639  int i;
1640 
1641 #ifndef VAPBSQUIET
1642  Vnm_tprint(1, "Destroying force arrays.\n");
1643 #endif
1644 
1645  for (i=0; i<nosh->ncalc; i++) {
1646 
1647  if (nforce[i] > 0) Vmem_free(mem, nforce[i], sizeof(AtomForce),
1648  (void **)&(atomForce[i]));
1649 
1650  }
1651 }
1652 
1653 VPUBLIC int writematMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg) {
1654 
1655  char writematstem[VMAX_ARGLEN];
1656  char outpath[VMAX_ARGLEN];
1657  char mxtype[3];
1658  int strlenmax;
1659 
1660  if (nosh->bogus) return 1;
1661 
1662 #ifdef HAVE_MPI_H
1663  strlenmax = VMAX_ARGLEN-14;
1664  if (strlen(pbeparm->writematstem) > strlenmax) {
1665  Vnm_tprint(2, " Matrix name (%s) too long (%d char max)!\n",
1666  pbeparm->writematstem, strlenmax);
1667  Vnm_tprint(2, " Not writing matrix!\n");
1668  return 0;
1669  }
1670  sprintf(writematstem, "%s-PE%d", pbeparm->writematstem, rank);
1671 #else
1672  strlenmax = (int)(VMAX_ARGLEN)-1;
1673  if ((int)strlen(pbeparm->writematstem) > strlenmax) {
1674  Vnm_tprint(2, " Matrix name (%s) too long (%d char max)!\n",
1675  pbeparm->writematstem, strlenmax);
1676  Vnm_tprint(2, " Not writing matrix!\n");
1677  return 0;
1678  }
1679  if(nosh->ispara == 1){
1680  sprintf(writematstem, "%s-PE%d", pbeparm->writematstem,nosh->proc_rank);
1681  }else{
1682  sprintf(writematstem, "%s", pbeparm->writematstem);
1683  }
1684 #endif
1685 
1686  if (pbeparm->writemat == 1) {
1687  strlenmax = VMAX_ARGLEN-5;
1688  if ((int)strlen(pbeparm->writematstem) > strlenmax) {
1689  Vnm_tprint(2, " Matrix name (%s) too long (%d char max)!\n",
1690  pbeparm->writematstem, strlenmax);
1691  Vnm_tprint(2, " Not writing matrix!\n");
1692  return 0;
1693  }
1694  sprintf(outpath, "%s.%s", writematstem, "mat");
1695  mxtype[0] = 'R';
1696  mxtype[1] = 'S';
1697  mxtype[2] = 'A';
1698  /* Poisson operator only */
1699  if (pbeparm->writematflag == 0) {
1700  Vnm_tprint( 1, " Writing Poisson operator matrix \
1701 to %s...\n", outpath);
1702 
1703  /* Linearization of Poisson-Boltzmann operator around solution */
1704  } else if (pbeparm->writematflag == 1) {
1705  Vnm_tprint( 1, " Writing linearization of full \
1706 Poisson-Boltzmann operator matrix to %s...\n", outpath);
1707 
1708  } else {
1709  Vnm_tprint( 2, " Bogus matrix specification\
1710 (%d)!\n", pbeparm->writematflag);
1711  return 0;
1712  }
1713 
1714  Vnm_tprint(0, " Printing operator...\n");
1715  //Vpmg_printColComp(pmg, outpath, outpath, mxtype,
1716  // pbeparm->writematflag);
1717  return 0;
1718 
1719  }
1720 
1721  return 1;
1722 }
1723 
1724 VPUBLIC void storeAtomEnergy(Vpmg *pmg, int icalc, double **atomEnergy,
1725  int *nenergy){
1726 
1727  Vatom *atom;
1728  Valist *alist;
1729  int i;
1730 
1731  alist = pmg->pbe->alist;
1732  *nenergy = Valist_getNumberAtoms(alist);
1733  *atomEnergy = (double *)Vmem_malloc(pmg->vmem, *nenergy, sizeof(double));
1734 
1735  for (i=0; i<*nenergy; i++) {
1736  atom = Valist_getAtom(alist, i);
1737  (*atomEnergy)[i] = Vpmg_qfAtomEnergy(pmg, atom);
1738  }
1739 }
1740 
1741 VPUBLIC int writedataFlat(
1742  NOsh *nosh,
1743  Vcom *com,
1744  const char *fname,
1745  double totEnergy[NOSH_MAXCALC],
1746  double qfEnergy[NOSH_MAXCALC],
1747  double qmEnergy[NOSH_MAXCALC],
1748  double dielEnergy[NOSH_MAXCALC],
1749  int nenergy[NOSH_MAXCALC],
1750  double *atomEnergy[NOSH_MAXCALC],
1751  int nforce[NOSH_MAXCALC],
1752  AtomForce *atomForce[NOSH_MAXCALC]) {
1753 
1754  FILE *file;
1755  time_t now;
1756  int ielec, icalc, i, j;
1757  char *timestring = VNULL;
1758  PBEparm *pbeparm = VNULL;
1759  MGparm *mgparm = VNULL;
1760  double conversion, ltenergy, gtenergy, scalar;
1761 
1762  if (nosh->bogus) return 1;
1763 
1764  /* Initialize some variables */
1765 
1766  icalc = 0;
1767 
1768  file = fopen(fname, "w");
1769  if (file == VNULL) {
1770  Vnm_print(2, "writedataFlat: Problem opening virtual socket %s\n",
1771  fname);
1772  return 0;
1773  }
1774 
1775  /* Strip the newline character from the date */
1776 
1777  now = time(VNULL);
1778  timestring = ctime(&now);
1779  fprintf(file,"%s\n", timestring);
1780 
1781  for (ielec=0; ielec<nosh->nelec;ielec++) { /* elec loop */
1782 
1783  /* Initialize per-elec pointers */
1784 
1785  mgparm = nosh->calc[icalc]->mgparm;
1786  pbeparm = nosh->calc[icalc]->pbeparm;
1787 
1788  /* Convert from kT/e to kJ/mol */
1789  conversion = Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na;
1790 
1791  fprintf(file,"elec");
1792  if (Vstring_strcasecmp(nosh->elecname[ielec], "") != 0) {
1793  fprintf(file," name %s\n", nosh->elecname[ielec]);
1794  } else fprintf(file, "\n");
1795 
1796  switch (mgparm->type) {
1797  case MCT_DUMMY:
1798  fprintf(file," mg-dummy\n");
1799  break;
1800  case MCT_MANUAL:
1801  fprintf(file," mg-manual\n");
1802  break;
1803  case MCT_AUTO:
1804  fprintf(file," mg-auto\n");
1805  break;
1806  case MCT_PARALLEL:
1807  fprintf(file," mg-para\n");
1808  break;
1809  default:
1810  break;
1811  }
1812 
1813  fprintf(file," mol %d\n", pbeparm->molid);
1814  fprintf(file," dime %d %d %d\n", mgparm->dime[0], mgparm->dime[1],\
1815  mgparm->dime[2]);
1816 
1817  switch (pbeparm->pbetype) {
1818  case PBE_NPBE:
1819  fprintf(file," npbe\n");
1820  break;
1821  case PBE_LPBE:
1822  fprintf(file," lpbe\n");
1823  break;
1824  default:
1825  break;
1826  }
1827 
1828  if (pbeparm->nion > 0) {
1829  for (i=0; i<pbeparm->nion; i++) {
1830  fprintf(file," ion %4.3f %4.3f %4.3f\n",
1831  pbeparm->ionr[i], pbeparm->ionq[i], pbeparm->ionc[i]);
1832  }
1833  }
1834 
1835  fprintf(file," pdie %4.3f\n", pbeparm->pdie);
1836  fprintf(file," sdie %4.3f\n", pbeparm->sdie);
1837 
1838  switch (pbeparm->srfm) {
1839  case 0:
1840  fprintf(file," srfm mol\n");
1841  fprintf(file," srad %4.3f\n", pbeparm->srad);
1842  break;
1843  case 1:
1844  fprintf(file," srfm smol\n");
1845  fprintf(file," srad %4.3f\n", pbeparm->srad);
1846  break;
1847  case 2:
1848  fprintf(file," srfm spl2\n");
1849  fprintf(file," srad %4.3f\n", pbeparm->srad);
1850  break;
1851  default:
1852  break;
1853  }
1854 
1855  switch (pbeparm->bcfl) {
1856  case BCFL_ZERO:
1857  fprintf(file," bcfl zero\n");
1858  break;
1859  case BCFL_SDH:
1860  fprintf(file," bcfl sdh\n");
1861  break;
1862  case BCFL_MDH:
1863  fprintf(file," bcfl mdh\n");
1864  break;
1865  case BCFL_FOCUS:
1866  fprintf(file," bcfl focus\n");
1867  break;
1868  case BCFL_MAP:
1869  fprintf(file," bcfl map\n");
1870  break;
1871  case BCFL_MEM:
1872  fprintf(file," bcfl mem\n");
1873  break;
1874  default:
1875  break;
1876  }
1877 
1878  fprintf(file," temp %4.3f\n", pbeparm->temp);
1879 
1880  for (;icalc<=nosh->elec2calc[ielec];icalc++){ /* calc loop */
1881 
1882  /* Reinitialize per-calc pointers */
1883  mgparm = nosh->calc[icalc]->mgparm;
1884  pbeparm = nosh->calc[icalc]->pbeparm;
1885 
1886  fprintf(file," calc\n");
1887  fprintf(file," id %i\n", (icalc+1));
1888  fprintf(file," grid %4.3f %4.3f %4.3f\n",
1889  mgparm->grid[0], mgparm->grid[1], mgparm->grid[2]);
1890  fprintf(file," glen %4.3f %4.3f %4.3f\n",
1891  mgparm->glen[0], mgparm->glen[1], mgparm->glen[2]);
1892 
1893  if (pbeparm->calcenergy == PCE_TOTAL) {
1894  fprintf(file," totEnergy %1.12E kJ/mol\n",
1895  (totEnergy[icalc]*conversion));
1896  } if (pbeparm->calcenergy == PCE_COMPS) {
1897  fprintf(file," totEnergy %1.12E kJ/mol\n",
1898  (totEnergy[icalc]*conversion));
1899  fprintf(file," qfEnergy %1.12E kJ/mol\n",
1900  (0.5*qfEnergy[icalc]*conversion));
1901  fprintf(file," qmEnergy %1.12E kJ/mol\n",
1902  (qmEnergy[icalc]*conversion));
1903  fprintf(file," dielEnergy %1.12E kJ/mol\n",
1904  (dielEnergy[icalc]*conversion));
1905  for (i=0; i<nenergy[icalc]; i++){
1906  fprintf(file," atom %i %1.12E kJ/mol\n", i,
1907  (0.5*atomEnergy[icalc][i]*conversion));
1908 
1909  }
1910  }
1911 
1912  if (pbeparm->calcforce == PCF_TOTAL) {
1913  fprintf(file," qfForce %1.12E %1.12E %1.12E kJ/mol/A\n",
1914  (atomForce[icalc][0].qfForce[0]*conversion),
1915  (atomForce[icalc][0].qfForce[1]*conversion),
1916  (atomForce[icalc][0].qfForce[2]*conversion));
1917  fprintf(file," ibForce %1.12E %1.12E %1.12E kJ/mol/A\n",
1918  (atomForce[icalc][0].ibForce[0]*conversion),
1919  (atomForce[icalc][0].ibForce[1]*conversion),
1920  (atomForce[icalc][0].ibForce[2]*conversion));
1921  fprintf(file," dbForce %1.12E %1.12E %1.12E kJ/mol/A\n",
1922  (atomForce[icalc][0].dbForce[0]*conversion),
1923  (atomForce[icalc][0].dbForce[1]*conversion),
1924  (atomForce[icalc][0].dbForce[2]*conversion));
1925  }
1926  fprintf(file," end\n");
1927  }
1928 
1929  fprintf(file,"end\n");
1930  }
1931 
1932 /* Handle print energy statements */
1933 
1934 for (i=0; i<nosh->nprint; i++) {
1935 
1936  if (nosh->printwhat[i] == NPT_ENERGY) {
1937 
1938  fprintf(file,"print energy");
1939  fprintf(file," %d", nosh->printcalc[i][0]+1);
1940 
1941  for (j=1; j<nosh->printnarg[i]; j++) {
1942  if (nosh->printop[i][j-1] == 0) fprintf(file," +");
1943  else if (nosh->printop[i][j-1] == 1) fprintf(file, " -");
1944  fprintf(file, " %d", nosh->printcalc[i][j]+1);
1945  }
1946 
1947  fprintf(file, "\n");
1948  icalc = nosh->elec2calc[nosh->printcalc[i][0]];
1949 
1950  ltenergy = Vunit_kb * (1e-3) * Vunit_Na * \
1951  nosh->calc[icalc]->pbeparm->temp * totEnergy[icalc];
1952 
1953  for (j=1; j<nosh->printnarg[i]; j++) {
1954  icalc = nosh->elec2calc[nosh->printcalc[i][j]];
1955  /* Add or subtract? */
1956  if (nosh->printop[i][j-1] == 0) scalar = 1.0;
1957  else if (nosh->printop[i][j-1] == 1) scalar = -1.0;
1958  /* Accumulate */
1959  ltenergy += (scalar * Vunit_kb * (1e-3) * Vunit_Na *
1960  nosh->calc[icalc]->pbeparm->temp * totEnergy[icalc]);
1961 
1962  Vcom_reduce(com, &ltenergy, &gtenergy, 1, 2, 0);
1963 
1964  }
1965  fprintf(file," localEnergy %1.12E kJ/mol\n", \
1966  ltenergy);
1967  fprintf(file," globalEnergy %1.12E kJ/mol\nend\n", \
1968  gtenergy);
1969  }
1970 }
1971 
1972 fclose(file);
1973 
1974 return 1;
1975 }
1976 
1977 VPUBLIC int writedataXML(NOsh *nosh, Vcom *com, const char *fname,
1978  double totEnergy[NOSH_MAXCALC],
1979  double qfEnergy[NOSH_MAXCALC],
1980  double qmEnergy[NOSH_MAXCALC],
1981  double dielEnergy[NOSH_MAXCALC],
1982  int nenergy[NOSH_MAXCALC],
1983  double *atomEnergy[NOSH_MAXCALC],
1984  int nforce[NOSH_MAXCALC],
1985  AtomForce *atomForce[NOSH_MAXCALC]) {
1986 
1987  FILE *file;
1988  time_t now;
1989  int ielec, icalc, i, j;
1990  char *timestring = VNULL;
1991  char *c = VNULL;
1992  PBEparm *pbeparm = VNULL;
1993  MGparm *mgparm = VNULL;
1994  double conversion, ltenergy, gtenergy, scalar;
1995 
1996  if (nosh->bogus) return 1;
1997 
1998  /* Initialize some variables */
1999 
2000  icalc = 0;
2001 
2002  file = fopen(fname, "w");
2003  if (file == VNULL) {
2004  Vnm_print(2, "writedataXML: Problem opening virtual socket %s\n",
2005  fname);
2006  return 0;
2007  }
2008 
2009  fprintf(file,"<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
2010  fprintf(file,"<APBS>\n");
2011 
2012  /* Strip the newline character from the date */
2013 
2014  now = time(VNULL);
2015  timestring = ctime(&now);
2016  for(c = timestring; *c != '\n'; c++);
2017  *c = '\0';
2018  fprintf(file," <date>%s</date>\n", timestring);
2019 
2020  for (ielec=0; ielec<nosh->nelec;ielec++){ /* elec loop */
2021 
2022  /* Initialize per-elec pointers */
2023 
2024  mgparm = nosh->calc[icalc]->mgparm;
2025  pbeparm = nosh->calc[icalc]->pbeparm;
2026 
2027  /* Convert from kT/e to kJ/mol */
2028  conversion = Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na;
2029 
2030  fprintf(file," <elec>\n");
2031  if (Vstring_strcasecmp(nosh->elecname[ielec], "") != 0) {
2032  fprintf(file," <name>%s</name>\n", nosh->elecname[ielec]);
2033  }
2034 
2035  switch (mgparm->type) {
2036  case MCT_DUMMY:
2037  fprintf(file," <type>mg-dummy</type>\n");
2038  break;
2039  case MCT_MANUAL:
2040  fprintf(file," <type>mg-manual</type>\n");
2041  break;
2042  case MCT_AUTO:
2043  fprintf(file," <type>mg-auto</type>\n");
2044  break;
2045  case MCT_PARALLEL:
2046  fprintf(file," <type>mg-para</type>\n");
2047  break;
2048  default:
2049  break;
2050  }
2051 
2052  fprintf(file," <molid>%d</molid>\n", pbeparm->molid);
2053  fprintf(file," <nx>%d</nx>\n", mgparm->dime[0]);
2054  fprintf(file," <ny>%d</ny>\n", mgparm->dime[1]);
2055  fprintf(file," <nz>%d</nz>\n", mgparm->dime[2]);
2056 
2057  switch (pbeparm->pbetype) {
2058  case PBE_NPBE:
2059  fprintf(file," <pbe>npbe</pbe>\n");
2060  break;
2061  case PBE_LPBE:
2062  fprintf(file," <pbe>lpbe</pbe>\n");
2063  break;
2064  default:
2065  break;
2066  }
2067 
2068  if (pbeparm->nion > 0) {
2069  for (i=0; i<pbeparm->nion; i++) {
2070  fprintf(file, " <ion>\n");
2071  fprintf(file," <radius>%4.3f A</radius>\n",
2072  pbeparm->ionr[i]);
2073  fprintf(file," <charge>%4.3f A</charge>\n",
2074  pbeparm->ionq[i]);
2075  fprintf(file," <concentration>%4.3f M</concentration>\n",
2076  pbeparm->ionc[i]);
2077  fprintf(file, " </ion>\n");
2078 
2079  }
2080  }
2081 
2082  fprintf(file," <pdie>%4.3f</pdie>\n", pbeparm->pdie);
2083  fprintf(file," <sdie>%4.3f</sdie>\n", pbeparm->sdie);
2084 
2085  switch (pbeparm->srfm) {
2086  case 0:
2087  fprintf(file," <srfm>mol</srfm>\n");
2088  fprintf(file," <srad>%4.3f</srad>\n", pbeparm->srad);
2089  break;
2090  case 1:
2091  fprintf(file," <srfm>smol</srfm>\n");
2092  fprintf(file," <srad>%4.3f</srad>\n", pbeparm->srad);
2093  break;
2094  case 2:
2095  fprintf(file," <srfm>spl2</srfm>\n");
2096  break;
2097  default:
2098  break;
2099  }
2100 
2101  switch (pbeparm->bcfl) {
2102  case BCFL_ZERO:
2103  fprintf(file," <bcfl>zero</bcfl>\n");
2104  break;
2105  case BCFL_SDH:
2106  fprintf(file," <bcfl>sdh</bcfl>\n");
2107  break;
2108  case BCFL_MDH:
2109  fprintf(file," <bcfl>mdh</bcfl>\n");
2110  break;
2111  case BCFL_FOCUS:
2112  fprintf(file," <bcfl>focus</bcfl>\n");
2113  break;
2114  case BCFL_MAP:
2115  fprintf(file," <bcfl>map</bcfl>\n");
2116  break;
2117  case BCFL_MEM:
2118  fprintf(file," <bcfl>mem</bcfl>\n");
2119  break;
2120  default:
2121  break;
2122  }
2123 
2124  fprintf(file," <temp>%4.3f K</temp>\n", pbeparm->temp);
2125 
2126  for (;icalc<=nosh->elec2calc[ielec];icalc++){ /* calc loop */
2127 
2128  /* Reinitialize per-calc pointers */
2129  mgparm = nosh->calc[icalc]->mgparm;
2130  pbeparm = nosh->calc[icalc]->pbeparm;
2131 
2132  fprintf(file," <calc>\n");
2133  fprintf(file," <id>%i</id>\n", (icalc+1));
2134  fprintf(file," <hx>%4.3f A</hx>\n", mgparm->grid[0]);
2135  fprintf(file," <hy>%4.3f A</hy>\n", mgparm->grid[1]);
2136  fprintf(file," <hz>%4.3f A</hz>\n", mgparm->grid[2]);
2137  fprintf(file," <xlen>%4.3f A</xlen>\n", mgparm->glen[0]);
2138  fprintf(file," <ylen>%4.3f A</ylen>\n", mgparm->glen[1]);
2139  fprintf(file," <zlen>%4.3f A</zlen>\n", mgparm->glen[2]);
2140 
2141  if (pbeparm->calcenergy == PCE_TOTAL) {
2142  fprintf(file," <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2143  (totEnergy[icalc]*conversion));
2144  } else if (pbeparm->calcenergy == PCE_COMPS) {
2145  fprintf(file," <totEnergy>%1.12E kJ/mol</totEnergy>\n",
2146  (totEnergy[icalc]*conversion));
2147  fprintf(file," <qfEnergy>%1.12E kJ/mol</qfEnergy>\n",
2148  (0.5*qfEnergy[icalc]*conversion));
2149  fprintf(file," <qmEnergy>%1.12E kJ/mol</qmEnergy>\n",
2150  (qmEnergy[icalc]*conversion));
2151  fprintf(file," <dielEnergy>%1.12E kJ/mol</dielEnergy>\n",
2152  (dielEnergy[icalc]*conversion));
2153  for (i=0; i<nenergy[icalc]; i++){
2154  fprintf(file," <atom>\n");
2155  fprintf(file," <id>%i</id>\n", i+1);
2156  fprintf(file," <energy>%1.12E kJ/mol</energy>\n",
2157  (0.5*atomEnergy[icalc][i]*conversion));
2158  fprintf(file," </atom>\n");
2159  }
2160  }
2161 
2162 
2163  if (pbeparm->calcforce == PCF_TOTAL) {
2164  fprintf(file," <qfforce_x>%1.12E</qfforce_x>\n",
2165  atomForce[icalc][0].qfForce[0]*conversion);
2166  fprintf(file," <qfforce_y>%1.12E</qfforce_y>\n",
2167  atomForce[icalc][0].qfForce[1]*conversion);
2168  fprintf(file," <qfforce_z>%1.12E</qfforce_z>\n",
2169  atomForce[icalc][0].qfForce[2]*conversion);
2170  fprintf(file," <ibforce_x>%1.12E</ibforce_x>\n",
2171  atomForce[icalc][0].ibForce[0]*conversion);
2172  fprintf(file," <ibforce_y>%1.12E</ibforce_y>\n",
2173  atomForce[icalc][0].ibForce[1]*conversion);
2174  fprintf(file," <ibforce_z>%1.12E</ibforce_z>\n",
2175  atomForce[icalc][0].ibForce[2]*conversion);
2176  fprintf(file," <dbforce_x>%1.12E</dbforce_x>\n",
2177  atomForce[icalc][0].dbForce[0]*conversion);
2178  fprintf(file," <dbforce_y>%1.12E</dbforce_y>\n",
2179  atomForce[icalc][0].dbForce[1]*conversion);
2180  fprintf(file," <dbforce_z>%1.12E</dbforce_z>\n",
2181  atomForce[icalc][0].dbForce[2]*conversion);
2182  }
2183 
2184  fprintf(file," </calc>\n");
2185  }
2186 
2187  fprintf(file," </elec>\n");
2188  }
2189 
2190 /* Handle print energy statements */
2191 
2192 for (i=0; i<nosh->nprint; i++) {
2193 
2194  if (nosh->printwhat[i] == NPT_ENERGY) {
2195 
2196  fprintf(file," <printEnergy>\n");
2197  fprintf(file," <equation>%d", nosh->printcalc[i][0]+1);
2198 
2199  for (j=1; j<nosh->printnarg[i]; j++) {
2200  if (nosh->printop[i][j-1] == 0) fprintf(file," +");
2201  else if (nosh->printop[i][j-1] == 1) fprintf(file, " -");
2202  fprintf(file, " %d", nosh->printcalc[i][j] +1);
2203  }
2204 
2205  fprintf(file, "</equation>\n");
2206  icalc = nosh->elec2calc[nosh->printcalc[i][0]];
2207 
2208  ltenergy = Vunit_kb * (1e-3) * Vunit_Na * \
2209  nosh->calc[icalc]->pbeparm->temp * totEnergy[icalc];
2210 
2211  for (j=1; j<nosh->printnarg[i]; j++) {
2212  icalc = nosh->elec2calc[nosh->printcalc[i][j]];
2213  /* Add or subtract? */
2214  if (nosh->printop[i][j-1] == 0) scalar = 1.0;
2215  else if (nosh->printop[i][j-1] == 1) scalar = -1.0;
2216  /* Accumulate */
2217  ltenergy += (scalar * Vunit_kb * (1e-3) * Vunit_Na *
2218  nosh->calc[icalc]->pbeparm->temp * totEnergy[icalc]);
2219  }
2220  Vcom_reduce(com, &ltenergy, &gtenergy, 1, 2, 0);
2221  fprintf(file," <localEnergy>%1.12E kJ/mol</localEnergy>\n", \
2222  ltenergy);
2223  fprintf(file," <globalEnergy>%1.12E kJ/mol</globalEnergy>\n", \
2224  gtenergy);
2225 
2226  fprintf(file," </printEnergy>\n");
2227  }
2228 }
2229 
2230 /* Add ending tags and close the file */
2231 fprintf(file,"</APBS>\n");
2232 fclose(file);
2233 
2234 return 1;
2235 }
2236 
2237 VPUBLIC int writedataMG(int rank,
2238  NOsh *nosh,
2239  PBEparm *pbeparm,
2240  Vpmg *pmg
2241  ) {
2242 
2243  char writestem[VMAX_ARGLEN];
2244  char outpath[VMAX_ARGLEN];
2245  char title[72];
2246  int i,
2247  nx,
2248  ny,
2249  nz,
2250  natoms;
2251  double hx,
2252  hy,
2253  hzed,
2254  xcent,
2255  ycent,
2256  zcent,
2257  xmin,
2258  ymin,
2259  zmin;
2260 
2261  Vgrid *grid;
2262  Vio *sock;
2263 
2264  if (nosh->bogus) return 1;
2265 
2266  for (i=0; i<pbeparm->numwrite; i++) {
2267 
2268  nx = pmg->pmgp->nx;
2269  ny = pmg->pmgp->ny;
2270  nz = pmg->pmgp->nz;
2271  hx = pmg->pmgp->hx;
2272  hy = pmg->pmgp->hy;
2273  hzed = pmg->pmgp->hzed;
2274 
2275  switch (pbeparm->writetype[i]) {
2276 
2277  case VDT_CHARGE:
2278 
2279  Vnm_tprint(1, " Writing charge distribution to ");
2280  xcent = pmg->pmgp->xcent;
2281  ycent = pmg->pmgp->ycent;
2282  zcent = pmg->pmgp->zcent;
2283  xmin = xcent - 0.5*(nx-1)*hx;
2284  ymin = ycent - 0.5*(ny-1)*hy;
2285  zmin = zcent - 0.5*(nz-1)*hzed;
2286  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_CHARGE, 0.0,
2287  pbeparm->pbetype, pbeparm));
2288  sprintf(title, "CHARGE DISTRIBUTION (e)");
2289  break;
2290 
2291  case VDT_POT:
2292 
2293  Vnm_tprint(1, " Writing potential to ");
2294  xcent = pmg->pmgp->xcent;
2295  ycent = pmg->pmgp->ycent;
2296  zcent = pmg->pmgp->zcent;
2297  xmin = xcent - 0.5*(nx-1)*hx;
2298  ymin = ycent - 0.5*(ny-1)*hy;
2299  zmin = zcent - 0.5*(nz-1)*hzed;
2300  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_POT, 0.0,
2301  pbeparm->pbetype, pbeparm));
2302  sprintf(title, "POTENTIAL (kT/e)");
2303  break;
2304 
2305  case VDT_SMOL:
2306 
2307  Vnm_tprint(1, " Writing molecular accessibility to ");
2308  xcent = pmg->pmgp->xcent;
2309  ycent = pmg->pmgp->ycent;
2310  zcent = pmg->pmgp->zcent;
2311  xmin = xcent - 0.5*(nx-1)*hx;
2312  ymin = ycent - 0.5*(ny-1)*hy;
2313  zmin = zcent - 0.5*(nz-1)*hzed;
2314  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_SMOL,
2315  pbeparm->srad, pbeparm->pbetype, pbeparm));
2316  sprintf(title,
2317  "SOLVENT ACCESSIBILITY -- MOLECULAR (%4.3f PROBE)",
2318  pbeparm->srad);
2319  break;
2320 
2321  case VDT_SSPL:
2322 
2323  Vnm_tprint(1, " Writing spline-based accessibility to ");
2324  xcent = pmg->pmgp->xcent;
2325  ycent = pmg->pmgp->ycent;
2326  zcent = pmg->pmgp->zcent;
2327  xmin = xcent - 0.5*(nx-1)*hx;
2328  ymin = ycent - 0.5*(ny-1)*hy;
2329  zmin = zcent - 0.5*(nz-1)*hzed;
2330  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_SSPL,
2331  pbeparm->swin, pbeparm->pbetype, pbeparm));
2332  sprintf(title,
2333  "SOLVENT ACCESSIBILITY -- SPLINE (%4.3f WINDOW)",
2334  pbeparm->swin);
2335  break;
2336 
2337  case VDT_VDW:
2338 
2339  Vnm_tprint(1, " Writing van der Waals accessibility to ");
2340  xcent = pmg->pmgp->xcent;
2341  ycent = pmg->pmgp->ycent;
2342  zcent = pmg->pmgp->zcent;
2343  xmin = xcent - 0.5*(nx-1)*hx;
2344  ymin = ycent - 0.5*(ny-1)*hy;
2345  zmin = zcent - 0.5*(nz-1)*hzed;
2346  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_VDW, 0.0,
2347  pbeparm->pbetype, pbeparm));
2348  sprintf(title, "SOLVENT ACCESSIBILITY -- VAN DER WAALS");
2349  break;
2350 
2351  case VDT_IVDW:
2352 
2353  Vnm_tprint(1, " Writing ion accessibility to ");
2354  xcent = pmg->pmgp->xcent;
2355  ycent = pmg->pmgp->ycent;
2356  zcent = pmg->pmgp->zcent;
2357  xmin = xcent - 0.5*(nx-1)*hx;
2358  ymin = ycent - 0.5*(ny-1)*hy;
2359  zmin = zcent - 0.5*(nz-1)*hzed;
2360  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_IVDW,
2361  pmg->pbe->maxIonRadius, pbeparm->pbetype, pbeparm));
2362  sprintf(title,
2363  "ION ACCESSIBILITY -- SPLINE (%4.3f RADIUS)",
2364  pmg->pbe->maxIonRadius);
2365  break;
2366 
2367  case VDT_LAP:
2368 
2369  Vnm_tprint(1, " Writing potential Laplacian to ");
2370  xcent = pmg->pmgp->xcent;
2371  ycent = pmg->pmgp->ycent;
2372  zcent = pmg->pmgp->zcent;
2373  xmin = xcent - 0.5*(nx-1)*hx;
2374  ymin = ycent - 0.5*(ny-1)*hy;
2375  zmin = zcent - 0.5*(nz-1)*hzed;
2376  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_LAP, 0.0,
2377  pbeparm->pbetype, pbeparm));
2378  sprintf(title,
2379  "POTENTIAL LAPLACIAN (kT/e/A^2)");
2380  break;
2381 
2382  case VDT_EDENS:
2383 
2384  Vnm_tprint(1, " Writing energy density to ");
2385  xcent = pmg->pmgp->xcent;
2386  ycent = pmg->pmgp->ycent;
2387  zcent = pmg->pmgp->zcent;
2388  xmin = xcent - 0.5*(nx-1)*hx;
2389  ymin = ycent - 0.5*(ny-1)*hy;
2390  zmin = zcent - 0.5*(nz-1)*hzed;
2391  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_EDENS, 0.0,
2392  pbeparm->pbetype, pbeparm));
2393  sprintf(title, "ENERGY DENSITY (kT/e/A)^2");
2394  break;
2395 
2396  case VDT_NDENS:
2397 
2398  Vnm_tprint(1, " Writing number density to ");
2399  xcent = pmg->pmgp->xcent;
2400  ycent = pmg->pmgp->ycent;
2401  zcent = pmg->pmgp->zcent;
2402  xmin = xcent - 0.5*(nx-1)*hx;
2403  ymin = ycent - 0.5*(ny-1)*hy;
2404  zmin = zcent - 0.5*(nz-1)*hzed;
2405  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_NDENS, 0.0,
2406  pbeparm->pbetype, pbeparm));
2407  sprintf(title,
2408  "ION NUMBER DENSITY (M)");
2409  break;
2410 
2411  case VDT_QDENS:
2412 
2413  Vnm_tprint(1, " Writing charge density to ");
2414  xcent = pmg->pmgp->xcent;
2415  ycent = pmg->pmgp->ycent;
2416  zcent = pmg->pmgp->zcent;
2417  xmin = xcent - 0.5*(nx-1)*hx;
2418  ymin = ycent - 0.5*(ny-1)*hy;
2419  zmin = zcent - 0.5*(nz-1)*hzed;
2420  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_QDENS, 0.0,
2421  pbeparm->pbetype, pbeparm));
2422  sprintf(title,
2423  "ION CHARGE DENSITY (e_c * M)");
2424  break;
2425 
2426  case VDT_DIELX:
2427 
2428  Vnm_tprint(1, " Writing x-shifted dielectric map to ");
2429  xcent = pmg->pmgp->xcent + 0.5*hx;
2430  ycent = pmg->pmgp->ycent;
2431  zcent = pmg->pmgp->zcent;
2432  xmin = xcent - 0.5*(nx-1)*hx;
2433  ymin = ycent - 0.5*(ny-1)*hy;
2434  zmin = zcent - 0.5*(nz-1)*hzed;
2435  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_DIELX, 0.0,
2436  pbeparm->pbetype, pbeparm));
2437  sprintf(title,
2438  "X-SHIFTED DIELECTRIC MAP");
2439  break;
2440 
2441  case VDT_DIELY:
2442 
2443  Vnm_tprint(1, " Writing y-shifted dielectric map to ");
2444  xcent = pmg->pmgp->xcent;
2445  ycent = pmg->pmgp->ycent + 0.5*hy;
2446  zcent = pmg->pmgp->zcent;
2447  xmin = xcent - 0.5*(nx-1)*hx;
2448  ymin = ycent - 0.5*(ny-1)*hy;
2449  zmin = zcent - 0.5*(nz-1)*hzed;
2450  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_DIELY, 0.0,
2451  pbeparm->pbetype, pbeparm));
2452  sprintf(title,
2453  "Y-SHIFTED DIELECTRIC MAP");
2454  break;
2455 
2456  case VDT_DIELZ:
2457 
2458  Vnm_tprint(1, " Writing z-shifted dielectric map to ");
2459  xcent = pmg->pmgp->xcent;
2460  ycent = pmg->pmgp->ycent;
2461  zcent = pmg->pmgp->zcent + 0.5*hzed;
2462  xmin = xcent - 0.5*(nx-1)*hx;
2463  ymin = ycent - 0.5*(ny-1)*hy;
2464  zmin = zcent - 0.5*(nz-1)*hzed;
2465  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_DIELZ, 0.0,
2466  pbeparm->pbetype, pbeparm));
2467  sprintf(title,
2468  "Z-SHIFTED DIELECTRIC MAP");
2469  break;
2470 
2471  case VDT_KAPPA:
2472 
2473  Vnm_tprint(1, " Writing kappa map to ");
2474  xcent = pmg->pmgp->xcent;
2475  ycent = pmg->pmgp->ycent;
2476  zcent = pmg->pmgp->zcent;
2477  xmin = xcent - 0.5*(nx-1)*hx;
2478  ymin = ycent - 0.5*(ny-1)*hy;
2479  zmin = zcent - 0.5*(nz-1)*hzed;
2480  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_KAPPA, 0.0,
2481  pbeparm->pbetype, pbeparm));
2482  sprintf(title,
2483  "KAPPA MAP");
2484  break;
2485 
2486  case VDT_ATOMPOT:
2487 
2488  Vnm_tprint(1, " Writing atom potentials to ");
2489  xcent = pmg->pmgp->xcent;
2490  ycent = pmg->pmgp->ycent;
2491  zcent = pmg->pmgp->zcent;
2492  xmin = xcent - 0.5*(nx-1)*hx;
2493  ymin = ycent - 0.5*(ny-1)*hy;
2494  zmin = zcent - 0.5*(nz-1)*hzed;
2495  VASSERT(Vpmg_fillArray(pmg, pmg->rwork, VDT_ATOMPOT, 0.0,
2496  pbeparm->pbetype, pbeparm));
2497  sprintf(title,
2498  "ATOM POTENTIALS");
2499  break;
2500  default:
2501 
2502  Vnm_tprint(2, "Invalid data type for writing!\n");
2503  return 0;
2504  }
2505 
2506 
2507 #ifdef HAVE_MPI_H
2508  sprintf(writestem, "%s-PE%d", pbeparm->writestem[i], rank);
2509 #else
2510  if(nosh->ispara){
2511  sprintf(writestem, "%s-PE%d", pbeparm->writestem[i],nosh->proc_rank);
2512  }else{
2513  sprintf(writestem, "%s", pbeparm->writestem[i]);
2514  }
2515 #endif
2516 
2517  switch (pbeparm->writefmt[i]) {
2518 
2519  case VDF_DX:
2520  sprintf(outpath, "%s.%s", writestem, "dx");
2521  Vnm_tprint(1, "%s\n", outpath);
2522  grid = Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2523  pmg->rwork);
2524  Vgrid_writeDX(grid, "FILE", "ASC", VNULL, outpath, title,
2525  pmg->pvec);
2526  Vgrid_dtor(&grid);
2527  break;
2528 
2529  case VDF_AVS:
2530  sprintf(outpath, "%s.%s", writestem, "ucd");
2531  Vnm_tprint(1, "%s\n", outpath);
2532  Vnm_tprint(2, "Sorry, AVS format isn't supported for \
2533 uniform meshes yet!\n");
2534  break;
2535 
2536  case VDF_MCSF:
2537  sprintf(outpath, "%s.%s", writestem, "mcsf");
2538  Vnm_tprint(1, "%s\n", outpath);
2539  Vnm_tprint(2, "Sorry, MCSF format isn't supported for \
2540  uniform meshes yet!\n");
2541  break;
2542 
2543  case VDF_UHBD:
2544  sprintf(outpath, "%s.%s", writestem, "grd");
2545  Vnm_tprint(1, "%s\n", outpath);
2546  grid = Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2547  pmg->rwork);
2548  Vgrid_writeUHBD(grid, "FILE", "ASC", VNULL, outpath, title,
2549  pmg->pvec);
2550  Vgrid_dtor(&grid);
2551  break;
2552 
2553  case VDF_GZ:
2554  sprintf(outpath, "%s.%s", writestem, "dx.gz");
2555  Vnm_tprint(1, "%s\n", outpath);
2556  grid = Vgrid_ctor(nx, ny, nz, hx, hy, hzed, xmin, ymin, zmin,
2557  pmg->rwork);
2558  Vgrid_writeGZ(grid, "FILE", "ASC", VNULL, outpath, title,
2559  pmg->pvec);
2560  Vgrid_dtor(&grid);
2561  break;
2562  case VDF_FLAT:
2563  sprintf(outpath, "%s.%s", writestem, "txt");
2564  Vnm_tprint(1, "%s\n", outpath);
2565  Vnm_print(0, "routines: Opening virtual socket...\n");
2566  sock = Vio_ctor("FILE","ASC",VNULL,outpath,"w");
2567  if (sock == VNULL) {
2568  Vnm_print(2, "routines: Problem opening virtual socket %s\n",
2569  outpath);
2570  return 0;
2571  }
2572  if (Vio_connect(sock, 0) < 0) {
2573  Vnm_print(2, "routines: Problem connecting virtual socket %s\n",
2574  outpath);
2575  return 0;
2576  }
2577  Vio_printf(sock, "# Data from %s\n", PACKAGE_STRING);
2578  Vio_printf(sock, "# \n");
2579  Vio_printf(sock, "# %s\n", title);
2580  Vio_printf(sock, "# \n");
2581  natoms = pmg->pbe->alist[pbeparm->molid-1].number;
2582  for(i=0;i<natoms;i++)
2583  Vio_printf(sock, "%12.6e\n", pmg->rwork[i]);
2584  break;
2585  default:
2586  Vnm_tprint(2, "Bogus data format (%d)!\n",
2587  pbeparm->writefmt[i]);
2588  break;
2589  }
2590 
2591  }
2592 
2593  return 1;
2594 }
2595 
2596 VPUBLIC double returnEnergy(Vcom *com,
2597  NOsh *nosh,
2598  double totEnergy[NOSH_MAXCALC],
2599  int iprint
2600  ){
2601 
2602  int iarg,
2603  calcid;
2604  double ltenergy,
2605  scalar;
2606 
2607  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
2608  if (nosh->calc[calcid]->pbeparm->calcenergy != PCE_NO) {
2609  ltenergy = Vunit_kb * (1e-3) * Vunit_Na *
2610  nosh->calc[calcid]->pbeparm->temp * totEnergy[calcid];
2611  } else {
2612  Vnm_tprint( 2, " No energy available in Calculation %d\n", calcid+1);
2613  return 0.0;
2614  }
2615  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++){
2616  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]];
2617  /* Add or substract */
2618  if (nosh->printop[iprint][iarg-1] == 0) scalar = 1.0;
2619  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
2620  /* Accumulate */
2621  ltenergy += (scalar * Vunit_kb * (1e-3) * Vunit_Na *
2622  nosh->calc[calcid]->pbeparm->temp * totEnergy[calcid]);
2623  }
2624 
2625  return ltenergy;
2626 }
2627 
2628 VPUBLIC int printEnergy(Vcom *com,
2629  NOsh *nosh,
2630  double totEnergy[NOSH_MAXCALC],
2631  int iprint
2632  ) {
2633 
2634  int iarg,
2635  calcid;
2636  double ltenergy,
2637  gtenergy,
2638  scalar;
2639 
2640  Vnm_tprint( 2, "Warning: The 'energy' print keyword is deprecated.\n" \
2641  " Use elecEnergy for electrostatics energy calcs.\n\n");
2642 
2643  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][0]], "") == 0){
2644  Vnm_tprint( 1, "print energy %d ", nosh->printcalc[iprint][0]+1);
2645  } else {
2646  Vnm_tprint( 1, "print energy %d (%s) ", nosh->printcalc[iprint][0]+1,
2647  nosh->elecname[nosh->printcalc[iprint][0]]);
2648  }
2649  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2650  if (nosh->printop[iprint][iarg-1] == 0)
2651  Vnm_tprint(1, "+ ");
2652  else if (nosh->printop[iprint][iarg-1] == 1)
2653  Vnm_tprint(1, "- ");
2654  else {
2655  Vnm_tprint( 2, "Undefined PRINT operation!\n");
2656  return 0;
2657  }
2658  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][iarg]],
2659  "") == 0) {
2660  Vnm_tprint( 1, "%d ", nosh->printcalc[iprint][iarg]+1);
2661  } else {
2662  Vnm_tprint( 1, "%d (%s) ", nosh->printcalc[iprint][iarg]+1,
2663  nosh->elecname[nosh->printcalc[iprint][iarg]]);
2664  }
2665  }
2666  Vnm_tprint(1, "end\n");
2667  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
2668  if (nosh->calc[calcid]->pbeparm->calcenergy != PCE_NO) {
2669  ltenergy = Vunit_kb * (1e-3) * Vunit_Na *
2670  nosh->calc[calcid]->pbeparm->temp * totEnergy[calcid];
2671  } else {
2672  Vnm_tprint( 2, " Didn't calculate energy in Calculation \
2673 #%d\n", calcid+1);
2674  return 0;
2675  }
2676  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2677  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]];
2678  /* Add or subtract? */
2679  if (nosh->printop[iprint][iarg-1] == 0) scalar = 1.0;
2680  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
2681  /* Accumulate */
2682  ltenergy += (scalar * Vunit_kb * (1e-3) * Vunit_Na *
2683  nosh->calc[calcid]->pbeparm->temp * totEnergy[calcid]);
2684  }
2685 
2686  Vnm_tprint( 1, " Local net energy (PE %d) = %1.12E kJ/mol\n",
2687  Vcom_rank(com), ltenergy);
2688  Vnm_tprint( 0, "printEnergy: Performing global reduction (sum)\n");
2689  Vcom_reduce(com, &ltenergy, &gtenergy, 1, 2, 0);
2690  Vnm_tprint( 1, " Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2691 
2692  return 1;
2693 
2694 }
2695 
2696 VPUBLIC int printElecEnergy(Vcom *com,
2697  NOsh *nosh,
2698  double totEnergy[NOSH_MAXCALC],
2699  int iprint
2700  ) {
2701 
2702  int iarg,
2703  calcid;
2704  double ltenergy,
2705  gtenergy,
2706  scalar;
2707 
2708  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][0]], "") == 0){
2709  Vnm_tprint( 1, "\nprint energy %d ", nosh->printcalc[iprint][0]+1);
2710  } else {
2711  Vnm_tprint( 1, "\nprint energy %d (%s) ", nosh->printcalc[iprint][0]+1,
2712  nosh->elecname[nosh->printcalc[iprint][0]]);
2713  }
2714  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2715  if (nosh->printop[iprint][iarg-1] == 0)
2716  Vnm_tprint(1, "+ ");
2717  else if (nosh->printop[iprint][iarg-1] == 1)
2718  Vnm_tprint(1, "- ");
2719  else {
2720  Vnm_tprint( 2, "Undefined PRINT operation!\n");
2721  return 0;
2722  }
2723  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][iarg]],
2724  "") == 0) {
2725  Vnm_tprint( 1, "%d ", nosh->printcalc[iprint][iarg]+1);
2726  } else {
2727  Vnm_tprint( 1, "%d (%s) ", nosh->printcalc[iprint][iarg]+1,
2728  nosh->elecname[nosh->printcalc[iprint][iarg]]);
2729  }
2730  }
2731  Vnm_tprint(1, "end\n");
2732  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
2733  if (nosh->calc[calcid]->pbeparm->calcenergy != PCE_NO) {
2734  ltenergy = Vunit_kb * (1e-3) * Vunit_Na *
2735  nosh->calc[calcid]->pbeparm->temp * totEnergy[calcid];
2736  } else {
2737  Vnm_tprint( 2, " Didn't calculate energy in Calculation \
2738 #%d\n", calcid+1);
2739  return 0;
2740  }
2741  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2742  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]];
2743  /* Add or subtract? */
2744  if (nosh->printop[iprint][iarg-1] == 0) scalar = 1.0;
2745  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
2746  /* Accumulate */
2747  ltenergy += (scalar * Vunit_kb * (1e-3) * Vunit_Na *
2748  nosh->calc[calcid]->pbeparm->temp * totEnergy[calcid]);
2749  }
2750 
2751  Vnm_tprint( 1, " Local net energy (PE %d) = %1.12E kJ/mol\n",
2752  Vcom_rank(com), ltenergy);
2753  Vnm_tprint( 0, "printEnergy: Performing global reduction (sum)\n");
2754  Vcom_reduce(com, &ltenergy, &gtenergy, 1, 2, 0);
2755  Vnm_tprint( 1, " Global net ELEC energy = %1.12E kJ/mol\n", gtenergy);
2756 
2757  return 1;
2758 
2759 }
2760 
2761 VPUBLIC int printApolEnergy(NOsh *nosh,
2762  int iprint
2763  ) {
2764 
2765  int iarg,
2766  calcid;
2767  double gtenergy,
2768  scalar;
2769  APOLparm *apolparm = VNULL;
2770 
2771  if (Vstring_strcasecmp(nosh->apolname[nosh->printcalc[iprint][0]], "") == 0){
2772  Vnm_tprint( 1, "\nprint APOL energy %d ", nosh->printcalc[iprint][0]+1);
2773  } else {
2774  Vnm_tprint( 1, "\nprint APOL energy %d (%s) ", nosh->printcalc[iprint][0]+1,
2775  nosh->apolname[nosh->printcalc[iprint][0]]);
2776  }
2777  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2778  if (nosh->printop[iprint][iarg-1] == 0)
2779  Vnm_tprint(1, "+ ");
2780  else if (nosh->printop[iprint][iarg-1] == 1)
2781  Vnm_tprint(1, "- ");
2782  else {
2783  Vnm_tprint( 2, "Undefined PRINT operation!\n");
2784  return 0;
2785  }
2786  if (Vstring_strcasecmp(nosh->apolname[nosh->printcalc[iprint][iarg]],
2787  "") == 0) {
2788  Vnm_tprint( 1, "%d ", nosh->printcalc[iprint][iarg]+1);
2789  } else {
2790  Vnm_tprint( 1, "%d (%s) ", nosh->printcalc[iprint][iarg]+1,
2791  nosh->apolname[nosh->printcalc[iprint][iarg]]);
2792  }
2793  }
2794  Vnm_tprint(1, "end\n");
2795 
2796  calcid = nosh->apol2calc[nosh->printcalc[iprint][0]];
2797  apolparm = nosh->calc[calcid]->apolparm;
2798 
2799  if (apolparm->calcenergy == ACE_TOTAL) {
2800  gtenergy = ((apolparm->gamma*apolparm->sasa) + (apolparm->press*apolparm->sav) + (apolparm->wcaEnergy));
2801  } else {
2802  Vnm_tprint( 2, " Didn't calculate energy in Calculation #%d\n", calcid+1);
2803  return 0;
2804  }
2805  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2806  calcid = nosh->apol2calc[nosh->printcalc[iprint][iarg]];
2807  apolparm = nosh->calc[calcid]->apolparm;
2808 
2809  /* Add or subtract? */
2810  if (nosh->printop[iprint][iarg-1] == 0) scalar = 1.0;
2811  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
2812  /* Accumulate */
2813  gtenergy += (scalar * ((apolparm->gamma*apolparm->sasa) +
2814  (apolparm->press*apolparm->sav) +
2815  (apolparm->wcaEnergy)));
2816 
2817  }
2818 
2819  Vnm_tprint( 1, " Global net APOL energy = %1.12E kJ/mol\n", gtenergy);
2820  return 1;
2821 }
2822 
2823 VPUBLIC int printForce(Vcom *com,
2824  NOsh *nosh,
2825  int nforce[NOSH_MAXCALC],
2826  AtomForce *atomForce[NOSH_MAXCALC],
2827  int iprint
2828  ) {
2829 
2830  int iarg,
2831  ifr,
2832  ivc,
2833  calcid,
2834  refnforce,
2835  refcalcforce;
2836  double temp,
2837  scalar,
2838  totforce[3];
2839  AtomForce *lforce,
2840  *gforce,
2841  *aforce;
2842 
2843  Vnm_tprint( 2, "Warning: The 'force' print keyword is deprecated.\n" \
2844  " Use elecForce for electrostatics force calcs.\n\n");
2845 
2846  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][0]], "") == 0){
2847  Vnm_tprint( 1, "print force %d ", nosh->printcalc[iprint][0]+1);
2848  } else {
2849  Vnm_tprint( 1, "print force %d (%s) ", nosh->printcalc[iprint][0]+1,
2850  nosh->elecname[nosh->printcalc[iprint][0]]);
2851  }
2852  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2853  if (nosh->printop[iprint][iarg-1] == 0)
2854  Vnm_tprint(1, "+ ");
2855  else if (nosh->printop[iprint][iarg-1] == 1)
2856  Vnm_tprint(1, "- ");
2857  else {
2858  Vnm_tprint( 2, "Undefined PRINT operation!\n");
2859  return 0;
2860  }
2861  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][iarg]],
2862  "") == 0) {
2863  Vnm_tprint( 1, "%d ", nosh->printcalc[iprint][iarg]+1);
2864  } else {
2865  Vnm_tprint( 1, "%d (%s) ", nosh->printcalc[iprint][iarg]+1,
2866  nosh->elecname[nosh->printcalc[iprint][iarg]]);
2867  }
2868  }
2869  Vnm_tprint(1, "end\n");
2870 
2871  /* First, go through and make sure we did the same type of force
2872  * evaluation in each of the requested calculations */
2873  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
2874  refnforce = nforce[calcid];
2875  refcalcforce = nosh->calc[calcid]->pbeparm->calcforce;
2876  if (refcalcforce == PCF_NO) {
2877  Vnm_tprint( 2, " Didn't calculate force in calculation \
2878 #%d\n", calcid+1);
2879  return 0;
2880  }
2881  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2882  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]-1];
2883  if (nosh->calc[calcid]->pbeparm->calcforce != refcalcforce) {
2884  Vnm_tprint(2, " Inconsistent calcforce declarations in \
2885 calculations %d and %d\n", nosh->elec2calc[nosh->printcalc[iprint][0]]+1,
2886  calcid+1);
2887  return 0;
2888  }
2889  if (nforce[calcid] != refnforce) {
2890  Vnm_tprint(2, " Inconsistent number of forces evaluated in \
2891 calculations %d and %d\n", nosh->elec2calc[nosh->printcalc[iprint][0]]+1,
2892  calcid+1);
2893  return 0;
2894  }
2895  }
2896 
2897  /* Now, allocate space to accumulate the forces */
2898  lforce = (AtomForce *)Vmem_malloc(VNULL, refnforce, sizeof(AtomForce));
2899  gforce = (AtomForce *)Vmem_malloc(VNULL, refnforce, sizeof(AtomForce));
2900 
2901  /* Now, accumulate the forces */
2902  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
2903  aforce = atomForce[calcid];
2904  temp = nosh->calc[calcid]->pbeparm->temp;
2905 
2906  /* Load up the first calculation */
2907  if (refcalcforce == PCF_TOTAL) {
2908  /* Set to total force */
2909  for (ivc=0; ivc<3; ivc++) {
2910  lforce[0].qfForce[ivc] =
2911  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].qfForce[ivc];
2912  lforce[0].ibForce[ivc] =
2913  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].ibForce[ivc];
2914  lforce[0].dbForce[ivc] =
2915  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].dbForce[ivc];
2916  }
2917  } else if (refcalcforce == PCF_COMPS) {
2918  for (ifr=0; ifr<refnforce; ifr++) {
2919  for (ivc=0; ivc<3; ivc++) {
2920  lforce[ifr].qfForce[ivc] =
2921  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].qfForce[ivc];
2922  lforce[ifr].ibForce[ivc] =
2923  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].ibForce[ivc];
2924  lforce[ifr].dbForce[ivc] =
2925  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].dbForce[ivc];
2926  }
2927  }
2928  }
2929 
2930  /* Load up the rest of the calculations */
2931  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
2932  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]];
2933  temp = nosh->calc[calcid]->pbeparm->temp;
2934  aforce = atomForce[calcid];
2935  /* Get operation */
2936  if (nosh->printop[iprint][iarg-1] == 0) scalar = +1.0;
2937  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
2938  else scalar = 0.0;
2939  /* Accumulate */
2940  if (refcalcforce == PCF_TOTAL) {
2941  /* Set to total force */
2942  for (ivc=0; ivc<3; ivc++) {
2943  lforce[0].qfForce[ivc] +=
2944  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].qfForce[ivc]);
2945  lforce[0].ibForce[ivc] +=
2946  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].ibForce[ivc]);
2947  lforce[0].dbForce[ivc] +=
2948  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].dbForce[ivc]);
2949  }
2950  } else if (refcalcforce == PCF_COMPS) {
2951  for (ifr=0; ifr<refnforce; ifr++) {
2952  for (ivc=0; ivc<3; ivc++) {
2953  lforce[ifr].qfForce[ivc] +=
2954  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].qfForce[ivc]);
2955  lforce[ifr].ibForce[ivc] +=
2956  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].ibForce[ivc]);
2957  lforce[ifr].dbForce[ivc] +=
2958  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].dbForce[ivc]);
2959  }
2960  }
2961  }
2962  }
2963 
2964  Vnm_tprint( 0, "printEnergy: Performing global reduction (sum)\n");
2965  for (ifr=0; ifr<refnforce; ifr++) {
2966  Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
2967  Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
2968  Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
2969  }
2970 
2971 #if 0
2972  if (refcalcforce == PCF_TOTAL) {
2973  Vnm_tprint( 1, " Local net fixed charge force = \
2974 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
2975  lforce[0].qfForce[1], lforce[0].qfForce[2]);
2976  Vnm_tprint( 1, " Local net ionic boundary force = \
2977 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
2978  lforce[0].ibForce[1], lforce[0].ibForce[2]);
2979  Vnm_tprint( 1, " Local net dielectric boundary force = \
2980 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
2981  lforce[0].dbForce[1], lforce[0].dbForce[2]);
2982  } else if (refcalcforce == PCF_COMPS) {
2983  for (ifr=0; ifr<refnforce; ifr++) {
2984  Vnm_tprint( 1, " Local fixed charge force \
2985 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
2986  lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
2987  Vnm_tprint( 1, " Local ionic boundary force \
2988 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
2989  lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
2990  Vnm_tprint( 1, " Local dielectric boundary force \
2991 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
2992  lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
2993  }
2994  }
2995 #endif
2996 
2997  if (refcalcforce == PCF_TOTAL) {
2998  Vnm_tprint( 1, " Printing net forces (kJ/mol/A).\n");
2999  Vnm_tprint( 1, " Legend:\n");
3000  Vnm_tprint( 1, " tot -- Total force\n");
3001  Vnm_tprint( 1, " qf -- Fixed charge force\n");
3002  Vnm_tprint( 1, " db -- Dielectric boundary force\n");
3003  Vnm_tprint( 1, " ib -- Ionic boundary force\n");
3004 
3005  for (ivc=0; ivc<3; ivc++) {
3006  totforce[ivc] =
3007  gforce[0].qfForce[ivc] + gforce[0].ibForce[ivc] \
3008  + gforce[0].dbForce[ivc];
3009  }
3010 
3011  Vnm_tprint( 1, " tot %1.12E %1.12E %1.12E\n", totforce[0],
3012  totforce[1], totforce[2]);
3013  Vnm_tprint( 1, " qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3014  gforce[0].qfForce[1], gforce[0].qfForce[2]);
3015  Vnm_tprint( 1, " ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3016  gforce[0].ibForce[1], gforce[0].ibForce[2]);
3017  Vnm_tprint( 1, " db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3018  gforce[0].dbForce[1], gforce[0].dbForce[2]);
3019 
3020  } else if (refcalcforce == PCF_COMPS) {
3021 
3022  Vnm_tprint( 1, " Printing per-atom forces (kJ/mol/A).\n");
3023  Vnm_tprint( 1, " Legend:\n");
3024  Vnm_tprint( 1, " tot n -- Total force for atom n\n");
3025  Vnm_tprint( 1, " qf n -- Fixed charge force for atom n\n");
3026  Vnm_tprint( 1, " db n -- Dielectric boundary force for atom n\n");
3027  Vnm_tprint( 1, " ib n -- Ionic boundary force for atom n\n");
3028  Vnm_tprint( 1, " tot all -- Total force for system\n");
3029 
3030  totforce[0] = 0.0;
3031  totforce[1] = 0.0;
3032  totforce[2] = 0.0;
3033 
3034  for (ifr=0; ifr<refnforce; ifr++) {
3035  Vnm_tprint( 1, " qf %d %1.12E %1.12E %1.12E\n", ifr,
3036  gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3037  gforce[ifr].qfForce[2]);
3038  Vnm_tprint( 1, " ib %d %1.12E %1.12E %1.12E\n", ifr,
3039  gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3040  gforce[ifr].ibForce[2]);
3041  Vnm_tprint( 1, " db %d %1.12E %1.12E %1.12E\n", ifr,
3042  gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3043  gforce[ifr].dbForce[2]);
3044  Vnm_tprint( 1, " tot %d %1.12E %1.12E %1.12E\n", ifr,
3045  (gforce[ifr].dbForce[0] \
3046  + gforce[ifr].ibForce[0] +
3047  gforce[ifr].qfForce[0]),
3048  (gforce[ifr].dbForce[1] \
3049  + gforce[ifr].ibForce[1] +
3050  gforce[ifr].qfForce[1]),
3051  (gforce[ifr].dbForce[2] \
3052  + gforce[ifr].ibForce[2] +
3053  gforce[ifr].qfForce[2]));
3054  for (ivc=0; ivc<3; ivc++) {
3055  totforce[ivc] += (gforce[ifr].dbForce[ivc] \
3056  + gforce[ifr].ibForce[ivc] \
3057  + gforce[ifr].qfForce[ivc]);
3058  }
3059  }
3060  Vnm_tprint( 1, " tot all %1.12E %1.12E %1.12E\n", totforce[0],
3061  totforce[1], totforce[2]);
3062  }
3063 
3064  Vmem_free(VNULL, refnforce, sizeof(AtomForce), (void **)(&lforce));
3065  Vmem_free(VNULL, refnforce, sizeof(AtomForce), (void **)(&gforce));
3066 
3067  return 1;
3068 
3069 }
3070 
3071 VPUBLIC int printElecForce(Vcom *com,
3072  NOsh *nosh,
3073  int nforce[NOSH_MAXCALC],
3074  AtomForce *atomForce[NOSH_MAXCALC],
3075  int iprint
3076  ) {
3077 
3078  int iarg,
3079  ifr,
3080  ivc,
3081  calcid,
3082  refnforce,
3083  refcalcforce;
3084  double temp,
3085  scalar,
3086  totforce[3];
3087  AtomForce *lforce, *gforce, *aforce;
3088 
3089  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][0]], "") == 0){
3090  Vnm_tprint( 1, "print force %d ", nosh->printcalc[iprint][0]+1);
3091  } else {
3092  Vnm_tprint( 1, "print force %d (%s) ", nosh->printcalc[iprint][0]+1,
3093  nosh->elecname[nosh->printcalc[iprint][0]]);
3094  }
3095  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
3096  if (nosh->printop[iprint][iarg-1] == 0)
3097  Vnm_tprint(1, "+ ");
3098  else if (nosh->printop[iprint][iarg-1] == 1)
3099  Vnm_tprint(1, "- ");
3100  else {
3101  Vnm_tprint( 2, "Undefined PRINT operation!\n");
3102  return 0;
3103  }
3104  if (Vstring_strcasecmp(nosh->elecname[nosh->printcalc[iprint][iarg]],
3105  "") == 0) {
3106  Vnm_tprint( 1, "%d ", nosh->printcalc[iprint][iarg]+1);
3107  } else {
3108  Vnm_tprint( 1, "%d (%s) ", nosh->printcalc[iprint][iarg]+1,
3109  nosh->elecname[nosh->printcalc[iprint][iarg]]);
3110  }
3111  }
3112  Vnm_tprint(1, "end\n");
3113 
3114  /* First, go through and make sure we did the same type of force
3115  * evaluation in each of the requested calculations */
3116  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
3117  refnforce = nforce[calcid];
3118  refcalcforce = nosh->calc[calcid]->pbeparm->calcforce;
3119  if (refcalcforce == PCF_NO) {
3120  Vnm_tprint( 2, " Didn't calculate force in calculation \
3121 #%d\n", calcid+1);
3122  return 0;
3123  }
3124  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
3125  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]-1];
3126  if (nosh->calc[calcid]->pbeparm->calcforce != refcalcforce) {
3127  Vnm_tprint(2, " Inconsistent calcforce declarations in \
3128 calculations %d and %d\n", nosh->elec2calc[nosh->printcalc[iprint][0]]+1,
3129  calcid+1);
3130  return 0;
3131  }
3132  if (nforce[calcid] != refnforce) {
3133  Vnm_tprint(2, " Inconsistent number of forces evaluated in \
3134 calculations %d and %d\n", nosh->elec2calc[nosh->printcalc[iprint][0]]+1,
3135  calcid+1);
3136  return 0;
3137  }
3138  }
3139 
3140  /* Now, allocate space to accumulate the forces */
3141  lforce = (AtomForce *)Vmem_malloc(VNULL, refnforce, sizeof(AtomForce));
3142  gforce = (AtomForce *)Vmem_malloc(VNULL, refnforce, sizeof(AtomForce));
3143 
3144  /* Now, accumulate the forces */
3145  calcid = nosh->elec2calc[nosh->printcalc[iprint][0]];
3146  aforce = atomForce[calcid];
3147  temp = nosh->calc[calcid]->pbeparm->temp;
3148 
3149  /* Load up the first calculation */
3150  if (refcalcforce == PCF_TOTAL) {
3151  /* Set to total force */
3152  for (ivc=0; ivc<3; ivc++) {
3153  lforce[0].qfForce[ivc] =
3154  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].qfForce[ivc];
3155  lforce[0].ibForce[ivc] =
3156  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].ibForce[ivc];
3157  lforce[0].dbForce[ivc] =
3158  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].dbForce[ivc];
3159  }
3160  } else if (refcalcforce == PCF_COMPS) {
3161  for (ifr=0; ifr<refnforce; ifr++) {
3162  for (ivc=0; ivc<3; ivc++) {
3163  lforce[ifr].qfForce[ivc] =
3164  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].qfForce[ivc];
3165  lforce[ifr].ibForce[ivc] =
3166  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].ibForce[ivc];
3167  lforce[ifr].dbForce[ivc] =
3168  Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].dbForce[ivc];
3169  }
3170  }
3171  }
3172 
3173  /* Load up the rest of the calculations */
3174  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
3175  calcid = nosh->elec2calc[nosh->printcalc[iprint][iarg]];
3176  temp = nosh->calc[calcid]->pbeparm->temp;
3177  aforce = atomForce[calcid];
3178  /* Get operation */
3179  if (nosh->printop[iprint][iarg-1] == 0) scalar = +1.0;
3180  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
3181  else scalar = 0.0;
3182  /* Accumulate */
3183  if (refcalcforce == PCF_TOTAL) {
3184  /* Set to total force */
3185  for (ivc=0; ivc<3; ivc++) {
3186  lforce[0].qfForce[ivc] +=
3187  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].qfForce[ivc]);
3188  lforce[0].ibForce[ivc] +=
3189  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].ibForce[ivc]);
3190  lforce[0].dbForce[ivc] +=
3191  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[0].dbForce[ivc]);
3192  }
3193  } else if (refcalcforce == PCF_COMPS) {
3194  for (ifr=0; ifr<refnforce; ifr++) {
3195  for (ivc=0; ivc<3; ivc++) {
3196  lforce[ifr].qfForce[ivc] +=
3197  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].qfForce[ivc]);
3198  lforce[ifr].ibForce[ivc] +=
3199  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].ibForce[ivc]);
3200  lforce[ifr].dbForce[ivc] +=
3201  (scalar*Vunit_kb*(1e-3)*Vunit_Na*temp*aforce[ifr].dbForce[ivc]);
3202  }
3203  }
3204  }
3205  }
3206 
3207  Vnm_tprint( 0, "printEnergy: Performing global reduction (sum)\n");
3208  for (ifr=0; ifr<refnforce; ifr++) {
3209  Vcom_reduce(com, lforce[ifr].qfForce, gforce[ifr].qfForce, 3, 2, 0);
3210  Vcom_reduce(com, lforce[ifr].ibForce, gforce[ifr].ibForce, 3, 2, 0);
3211  Vcom_reduce(com, lforce[ifr].dbForce, gforce[ifr].dbForce, 3, 2, 0);
3212  }
3213 
3214 #if 0
3215  if (refcalcforce == PCF_TOTAL) {
3216  Vnm_tprint( 1, " Local net fixed charge force = \
3217 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].qfForce[0],
3218  lforce[0].qfForce[1], lforce[0].qfForce[2]);
3219  Vnm_tprint( 1, " Local net ionic boundary force = \
3220 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].ibForce[0],
3221  lforce[0].ibForce[1], lforce[0].ibForce[2]);
3222  Vnm_tprint( 1, " Local net dielectric boundary force = \
3223 (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", lforce[0].dbForce[0],
3224  lforce[0].dbForce[1], lforce[0].dbForce[2]);
3225  } else if (refcalcforce == PCF_COMPS) {
3226  for (ifr=0; ifr<refnforce; ifr++) {
3227  Vnm_tprint( 1, " Local fixed charge force \
3228 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].qfForce[0],
3229  lforce[ifr].qfForce[1], lforce[ifr].qfForce[2]);
3230  Vnm_tprint( 1, " Local ionic boundary force \
3231 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].ibForce[0],
3232  lforce[ifr].ibForce[1], lforce[ifr].ibForce[2]);
3233  Vnm_tprint( 1, " Local dielectric boundary force \
3234 (atom %d) = (%1.12E, %1.12E, %1.12E) kJ/mol/A\n", ifr, lforce[ifr].dbForce[0],
3235  lforce[ifr].dbForce[1], lforce[ifr].dbForce[2]);
3236  }
3237  }
3238 #endif
3239 
3240  if (refcalcforce == PCF_TOTAL) {
3241  Vnm_tprint( 1, " Printing net forces (kJ/mol/A).\n");
3242  Vnm_tprint( 1, " Legend:\n");
3243  Vnm_tprint( 1, " tot -- Total force\n");
3244  Vnm_tprint( 1, " qf -- Fixed charge force\n");
3245  Vnm_tprint( 1, " db -- Dielectric boundary force\n");
3246  Vnm_tprint( 1, " ib -- Ionic boundary force\n");
3247 
3248  for (ivc=0; ivc<3; ivc++) {
3249  totforce[ivc] =
3250  gforce[0].qfForce[ivc] + gforce[0].ibForce[ivc] \
3251  + gforce[0].dbForce[ivc];
3252  }
3253 
3254  Vnm_tprint( 1, " tot %1.12E %1.12E %1.12E\n", totforce[0],
3255  totforce[1], totforce[2]);
3256  Vnm_tprint( 1, " qf %1.12E %1.12E %1.12E\n", gforce[0].qfForce[0],
3257  gforce[0].qfForce[1], gforce[0].qfForce[2]);
3258  Vnm_tprint( 1, " ib %1.12E %1.12E %1.12E\n", gforce[0].ibForce[0],
3259  gforce[0].ibForce[1], gforce[0].ibForce[2]);
3260  Vnm_tprint( 1, " db %1.12E %1.12E %1.12E\n", gforce[0].dbForce[0],
3261  gforce[0].dbForce[1], gforce[0].dbForce[2]);
3262 
3263  } else if (refcalcforce == PCF_COMPS) {
3264 
3265  Vnm_tprint( 1, " Printing per-atom forces (kJ/mol/A).\n");
3266  Vnm_tprint( 1, " Legend:\n");
3267  Vnm_tprint( 1, " tot n -- Total force for atom n\n");
3268  Vnm_tprint( 1, " qf n -- Fixed charge force for atom n\n");
3269  Vnm_tprint( 1, " db n -- Dielectric boundary force for atom n\n");
3270  Vnm_tprint( 1, " ib n -- Ionic boundary force for atom n\n");
3271  Vnm_tprint( 1, " tot all -- Total force for system\n");
3272 
3273  totforce[0] = 0.0;
3274  totforce[1] = 0.0;
3275  totforce[2] = 0.0;
3276 
3277  for (ifr=0; ifr<refnforce; ifr++) {
3278  Vnm_tprint( 1, " qf %d %1.12E %1.12E %1.12E\n", ifr,
3279  gforce[ifr].qfForce[0], gforce[ifr].qfForce[1],
3280  gforce[ifr].qfForce[2]);
3281  Vnm_tprint( 1, " ib %d %1.12E %1.12E %1.12E\n", ifr,
3282  gforce[ifr].ibForce[0], gforce[ifr].ibForce[1],
3283  gforce[ifr].ibForce[2]);
3284  Vnm_tprint( 1, " db %d %1.12E %1.12E %1.12E\n", ifr,
3285  gforce[ifr].dbForce[0], gforce[ifr].dbForce[1],
3286  gforce[ifr].dbForce[2]);
3287  Vnm_tprint( 1, " tot %d %1.12E %1.12E %1.12E\n", ifr,
3288  (gforce[ifr].dbForce[0] \
3289  + gforce[ifr].ibForce[0] +
3290  gforce[ifr].qfForce[0]),
3291  (gforce[ifr].dbForce[1] \
3292  + gforce[ifr].ibForce[1] +
3293  gforce[ifr].qfForce[1]),
3294  (gforce[ifr].dbForce[2] \
3295  + gforce[ifr].ibForce[2] +
3296  gforce[ifr].qfForce[2]));
3297  for (ivc=0; ivc<3; ivc++) {
3298  totforce[ivc] += (gforce[ifr].dbForce[ivc] \
3299  + gforce[ifr].ibForce[ivc] \
3300  + gforce[ifr].qfForce[ivc]);
3301  }
3302  }
3303  Vnm_tprint( 1, " tot all %1.12E %1.12E %1.12E\n", totforce[0],
3304  totforce[1], totforce[2]);
3305  }
3306 
3307  Vmem_free(VNULL, refnforce, sizeof(AtomForce), (void **)(&lforce));
3308  Vmem_free(VNULL, refnforce, sizeof(AtomForce), (void **)(&gforce));
3309 
3310  return 1;
3311 
3312 }
3313 
3314 VPUBLIC int printApolForce(Vcom *com,
3315  NOsh *nosh,
3316  int nforce[NOSH_MAXCALC],
3317  AtomForce *atomForce[NOSH_MAXCALC],
3318  int iprint
3319  ) {
3320 
3321  int iarg,
3322  ifr,
3323  ivc,
3324  calcid,
3325  refnforce,
3326  refcalcforce;
3327  double temp,
3328  scalar,
3329  totforce[3];
3330  AtomForce *lforce,
3331  *gforce,
3332  *aforce;
3333 
3334  if (Vstring_strcasecmp(nosh->apolname[nosh->printcalc[iprint][0]], "") == 0){
3335  Vnm_tprint( 1, "\nprint APOL force %d ", nosh->printcalc[iprint][0]+1);
3336  } else {
3337  Vnm_tprint( 1, "\nprint APOL force %d (%s) ", nosh->printcalc[iprint][0]+1,
3338  nosh->apolname[nosh->printcalc[iprint][0]]);
3339  }
3340  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
3341  if (nosh->printop[iprint][iarg-1] == 0)
3342  Vnm_tprint(1, "+ ");
3343  else if (nosh->printop[iprint][iarg-1] == 1)
3344  Vnm_tprint(1, "- ");
3345  else {
3346  Vnm_tprint( 2, "Undefined PRINT operation!\n");
3347  return 0;
3348  }
3349  if (Vstring_strcasecmp(nosh->apolname[nosh->printcalc[iprint][iarg]],
3350  "") == 0) {
3351  Vnm_tprint( 1, "%d ", nosh->printcalc[iprint][iarg]+1);
3352  } else {
3353  Vnm_tprint( 1, "%d (%s) ", nosh->printcalc[iprint][iarg]+1,
3354  nosh->apolname[nosh->printcalc[iprint][iarg]]);
3355  }
3356  }
3357  Vnm_tprint(1, "end\n");
3358 
3359  /* First, go through and make sure we did the same type of force
3360  * evaluation in each of the requested calculations */
3361  calcid = nosh->apol2calc[nosh->printcalc[iprint][0]];
3362  refnforce = nforce[calcid];
3363  refcalcforce = nosh->calc[calcid]->apolparm->calcforce;
3364  if (refcalcforce == ACF_NO) {
3365  Vnm_tprint( 2, " Didn't calculate force in calculation \
3366 #%d\n", calcid+1);
3367  return 0;
3368  }
3369  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
3370  calcid = nosh->apol2calc[nosh->printcalc[iprint][iarg]-1];
3371  if (nosh->calc[calcid]->apolparm->calcforce != refcalcforce) {
3372  Vnm_tprint(2, " Inconsistent calcforce declarations in \
3373 calculations %d and %d\n", nosh->apol2calc[nosh->printcalc[iprint][0]]+1,
3374  calcid+1);
3375  return 0;
3376  }
3377  if (nforce[calcid] != refnforce) {
3378  Vnm_tprint(2, " Inconsistent number of forces evaluated in \
3379 calculations %d and %d\n", nosh->apol2calc[nosh->printcalc[iprint][0]]+1,
3380  calcid+1);
3381  return 0;
3382  }
3383  }
3384 
3385  /* Now, allocate space to accumulate the forces */
3386  lforce = (AtomForce *)Vmem_malloc(VNULL, refnforce, sizeof(AtomForce));
3387  gforce = (AtomForce *)Vmem_malloc(VNULL, refnforce, sizeof(AtomForce));
3388 
3389  /* Now, accumulate the forces */
3390  calcid = nosh->apol2calc[nosh->printcalc[iprint][0]];
3391  aforce = atomForce[calcid];
3392  temp = nosh->calc[calcid]->apolparm->temp;
3393 
3394  /* Load up the first calculation */
3395  if (refcalcforce == ACF_TOTAL) {
3396  /* Set to total force */
3397  for (ivc=0; ivc<3; ivc++) {
3398  lforce[0].sasaForce[ivc] = aforce[0].sasaForce[ivc];
3399  lforce[0].savForce[ivc] = aforce[0].savForce[ivc];
3400  lforce[0].wcaForce[ivc] = aforce[0].wcaForce[ivc];
3401  }
3402  } else if (refcalcforce == ACF_COMPS) {
3403  for (ifr=0; ifr<refnforce; ifr++) {
3404  for (ivc=0; ivc<3; ivc++) {
3405  lforce[ifr].sasaForce[ivc] = aforce[ifr].sasaForce[ivc];
3406  lforce[ifr].savForce[ivc] = aforce[ifr].savForce[ivc];
3407  lforce[ifr].wcaForce[ivc] = aforce[ifr].wcaForce[ivc];
3408  }
3409  }
3410  }
3411 
3412  /* Load up the rest of the calculations */
3413  for (iarg=1; iarg<nosh->printnarg[iprint]; iarg++) {
3414  calcid = nosh->apol2calc[nosh->printcalc[iprint][iarg]];
3415  temp = nosh->calc[calcid]->apolparm->temp;
3416  aforce = atomForce[calcid];
3417  /* Get operation */
3418  if (nosh->printop[iprint][iarg-1] == 0) scalar = +1.0;
3419  else if (nosh->printop[iprint][iarg-1] == 1) scalar = -1.0;
3420  else scalar = 0.0;
3421  /* Accumulate */
3422  if (refcalcforce == ACF_TOTAL) {
3423  /* Set to total force */
3424  for (ivc=0; ivc<3; ivc++) {
3425  lforce[0].sasaForce[ivc] += aforce[0].sasaForce[ivc];
3426  lforce[0].savForce[ivc] += aforce[0].savForce[ivc];
3427  lforce[0].wcaForce[ivc] += aforce[0].wcaForce[ivc];
3428  }
3429  } else if (refcalcforce == ACF_COMPS) {
3430  for (ifr=0; ifr<refnforce; ifr++) {
3431  for (ivc=0; ivc<3; ivc++) {
3432  lforce[ifr].sasaForce[ivc] += aforce[ifr].sasaForce[ivc];
3433  lforce[ifr].savForce[ivc] += aforce[ifr].savForce[ivc];
3434  lforce[ifr].wcaForce[ivc] += aforce[ifr].wcaForce[ivc];
3435  }
3436  }
3437  }
3438  }
3439 
3440  Vnm_tprint( 0, "printForce: Performing global reduction (sum)\n");
3441  for (ifr=0; ifr<refnforce; ifr++) {
3442  Vcom_reduce(com, lforce[ifr].sasaForce, gforce[ifr].sasaForce, 3, 2, 0);
3443  Vcom_reduce(com, lforce[ifr].savForce, gforce[ifr].savForce, 3, 2, 0);
3444  Vcom_reduce(com, lforce[ifr].wcaForce, gforce[ifr].wcaForce, 3, 2, 0);
3445  }
3446 
3447  if (refcalcforce == ACF_TOTAL) {
3448  Vnm_tprint( 1, " Printing net forces (kJ/mol/A)\n");
3449  Vnm_tprint( 1, " Legend:\n");
3450  Vnm_tprint( 1, " tot -- Total force\n");
3451  Vnm_tprint( 1, " sasa -- SASA force\n");
3452  Vnm_tprint( 1, " sav -- SAV force\n");
3453  Vnm_tprint( 1, " wca -- WCA force\n\n");
3454 
3455  for (ivc=0; ivc<3; ivc++) {
3456  totforce[ivc] =
3457  gforce[0].sasaForce[ivc] + gforce[0].savForce[ivc] \
3458  + gforce[0].wcaForce[ivc];
3459  }
3460 
3461  Vnm_tprint( 1, " tot %1.12E %1.12E %1.12E\n", totforce[0],
3462  totforce[1], totforce[2]);
3463  Vnm_tprint( 1, " sasa %1.12E %1.12E %1.12E\n", gforce[0].sasaForce[0],
3464  gforce[0].sasaForce[1], gforce[0].sasaForce[2]);
3465  Vnm_tprint( 1, " sav %1.12E %1.12E %1.12E\n", gforce[0].savForce[0],
3466  gforce[0].savForce[1], gforce[0].savForce[2]);
3467  Vnm_tprint( 1, " wca %1.12E %1.12E %1.12E\n", gforce[0].wcaForce[0],
3468  gforce[0].wcaForce[1], gforce[0].wcaForce[2]);
3469 
3470  } else if (refcalcforce == ACF_COMPS) {
3471 
3472  Vnm_tprint( 1, " Printing per atom forces (kJ/mol/A)\n");
3473  Vnm_tprint( 1, " Legend:\n");
3474  Vnm_tprint( 1, " tot n -- Total force for atom n\n");
3475  Vnm_tprint( 1, " sasa n -- SASA force for atom n\n");
3476  Vnm_tprint( 1, " sav n -- SAV force for atom n\n");
3477  Vnm_tprint( 1, " wca n -- WCA force for atom n\n");
3478  Vnm_tprint( 1, " tot all -- Total force for system\n");
3479 
3480  //Vnm_tprint( 1, " gamma, pressure, bconc are: %f %f %f\n\n",
3481  // gamma,press,bconc);
3482 
3483  totforce[0] = 0.0;
3484  totforce[1] = 0.0;
3485  totforce[2] = 0.0;
3486 
3487  for (ifr=0; ifr<refnforce; ifr++) {
3488  Vnm_tprint( 1, " sasa %d %1.12E %1.12E %1.12E\n", ifr,
3489  gforce[ifr].sasaForce[0], gforce[ifr].sasaForce[1],
3490  gforce[ifr].sasaForce[2]);
3491  Vnm_tprint( 1, " sav %d %1.12E %1.12E %1.12E\n", ifr,
3492  gforce[ifr].savForce[0], gforce[ifr].savForce[1],
3493  gforce[ifr].savForce[2]);
3494  Vnm_tprint( 1, " wca %d %1.12E %1.12E %1.12E\n", ifr,
3495  gforce[ifr].wcaForce[0], gforce[ifr].wcaForce[1],
3496  gforce[ifr].wcaForce[2]);
3497  Vnm_tprint( 1, " tot %d %1.12E %1.12E %1.12E\n", ifr,
3498  (gforce[ifr].wcaForce[0] \
3499  + gforce[ifr].savForce[0] +
3500  gforce[ifr].sasaForce[0]),
3501  (gforce[ifr].wcaForce[1] \
3502  + gforce[ifr].savForce[1] +
3503  gforce[ifr].sasaForce[1]),
3504  (gforce[ifr].wcaForce[2] \
3505  + gforce[ifr].savForce[2] +
3506  gforce[ifr].sasaForce[2]));
3507  for (ivc=0; ivc<3; ivc++) {
3508  totforce[ivc] += (gforce[ifr].wcaForce[ivc] \
3509  + gforce[ifr].savForce[ivc] \
3510  + gforce[ifr].sasaForce[ivc]);
3511  }
3512  }
3513  Vnm_tprint( 1, " tot all %1.12E %1.12E %1.12E\n", totforce[0],
3514  totforce[1], totforce[2]);
3515  }
3516 
3517  Vmem_free(VNULL, refnforce, sizeof(AtomForce), (void **)(&lforce));
3518  Vmem_free(VNULL, refnforce, sizeof(AtomForce), (void **)(&gforce));
3519 
3520  return 1;
3521 }
3522 
3523 #ifdef HAVE_MC_H
3524 
3525 
3526 VPUBLIC void killFE(NOsh *nosh,
3527  Vpbe *pbe[NOSH_MAXCALC],
3528  Vfetk *fetk[NOSH_MAXCALC],
3529  Gem *gm[NOSH_MAXMOL]
3530  ) {
3531 
3532  int i;
3533 
3534 #ifndef VAPBSQUIET
3535  Vnm_tprint(1, "Destroying finite element structures.\n");
3536 #endif
3537 
3538  for(i=0;i<nosh->ncalc;i++){
3539  Vpbe_dtor(&(pbe[i]));
3540  Vfetk_dtor(&(fetk[i]));
3541  }
3542  for (i=0; i<nosh->nmesh; i++) {
3543  Gem_dtor(&(gm[i]));
3544 }
3545 }
3546 
3554 VPUBLIC Vrc_Codes initFE(int icalc,
3555  NOsh *nosh,
3556  FEMparm *feparm,
3557  PBEparm *pbeparm,
3558  Vpbe *pbe[NOSH_MAXCALC],
3559  Valist *alist[NOSH_MAXMOL],
3560  Vfetk *fetk[NOSH_MAXCALC]
3561  ) {
3562 
3563  int iatom,
3564  imesh,
3565  i,
3566  j,
3567  theMol,
3568  focusFlag = 0;
3569  Vio *sock = VNULL;
3570  size_t bytesTotal,
3571  highWater;
3572  Vfetk_MeshLoad meshType;
3573  double length[3],
3574  center[3],
3575  sparm,
3576  q,
3577  iparm = 0.0;
3578  Vrc_Codes vrc;
3579  Valist *myalist;
3580  Vatom *atom = VNULL;
3582  Vnm_tstart(27, "Setup timer");
3583 
3584  /* Print some warning messages */
3585  if (pbeparm->useDielMap) Vnm_tprint(2, "FEM ignoring dielectric map!\n");
3586  if (pbeparm->useKappaMap) Vnm_tprint(2, "FEM ignoring kappa map!\n");
3587  if (pbeparm->useChargeMap) Vnm_tprint(2, "FEM ignoring charge map!\n");
3588 
3589  /* Fix mesh center for "GCENT MOL #" types of declarations. */
3590  Vnm_tprint(0, "Re-centering mesh...\n");
3591  theMol = pbeparm->molid-1;
3592  myalist = alist[theMol];
3593  for (j=0; j<3; j++) {
3594  if (theMol < nosh->nmol) {
3595  center[j] = (myalist)->center[j];
3596  } else{
3597  Vnm_tprint(2, "ERROR! Bogus molecule number (%d)!\n",
3598  (theMol+1));
3599  return VRC_FAILURE;
3600  }
3601  }
3602 
3603  /* Check for completely-neutral molecule */
3604  q = 0;
3605  for (iatom=0; iatom<Valist_getNumberAtoms(myalist); iatom++) {
3606  atom = Valist_getAtom(myalist, iatom);
3607  q += VSQR(Vatom_getCharge(atom));
3608  }
3609  /* D. Gohara 10/22/09 - disabled
3610  if (q < (1e-6)) {
3611  Vnm_tprint(2, "Molecule #%d is uncharged!\n", pbeparm->molid);
3612  Vnm_tprint(2, "Sum square charge = %g!\n", q);
3613  return VRC_FAILURE;
3614  }
3615  */
3616 
3617  /* Set the femparm pkey value based on the presence of an HB solver */
3618 #ifdef USE_HB
3619  feparm->pkey = 1;
3620 #else
3621  feparm->pkey = 0;
3622 #endif
3623 
3624  /* Set up PBE object */
3625  Vnm_tprint(0, "Setting up PBE object...\n");
3626  if (pbeparm->srfm == VSM_SPLINE) {
3627  sparm = pbeparm->swin;
3628  }
3629  else {
3630  sparm = pbeparm->srad;
3631  }
3632  if (pbeparm->nion > 0) {
3633  iparm = pbeparm->ionr[0];
3634  }
3635 
3636  pbe[icalc] = Vpbe_ctor(myalist, pbeparm->nion,
3637  pbeparm->ionc, pbeparm->ionr, pbeparm->ionq,
3638  pbeparm->temp, pbeparm->pdie,
3639  pbeparm->sdie, sparm, focusFlag, pbeparm->sdens,
3640  pbeparm->zmem, pbeparm->Lmem, pbeparm->mdie,
3641  pbeparm->memv);
3642 
3643  /* Print a few derived parameters */
3644  Vnm_tprint(1, " Debye length: %g A\n", Vpbe_getDeblen(pbe[icalc]));
3645 
3646  /* Set up FEtk objects */
3647  Vnm_tprint(0, "Setting up FEtk object...\n");
3648  fetk[icalc] = Vfetk_ctor(pbe[icalc], pbeparm->pbetype);
3649  Vfetk_setParameters(fetk[icalc], pbeparm, feparm);
3650 
3651  /* Build mesh - this merely loads an MCSF file from an external source if one is specified or uses the
3652  * current molecule and sets center/length values based on that molecule if no external source is
3653  * specified. */
3654  Vnm_tprint(0, "Setting up mesh...\n");
3655  sock = VNULL;
3656  if (feparm->useMesh) {
3657  imesh = feparm->meshID-1;
3658  meshType = VML_EXTERNAL;
3659  for (i=0; i<3; i++) {
3660  center[i] = 0.0;
3661  length[i] = 0.0;
3662  }
3663  Vnm_print(0, "Using mesh %d (%s) in calculation.\n", imesh+1,
3664  nosh->meshpath[imesh]);
3665  switch (nosh->meshfmt[imesh]) {
3666  case VDF_DX:
3667  Vnm_tprint(2, "DX finite element mesh input not supported yet!\n");
3668  return VRC_FAILURE;
3669  case VDF_UHBD:
3670  Vnm_tprint( 2, "UHBD finite element mesh input not supported!\n");
3671  return VRC_FAILURE;
3672  case VDF_AVS:
3673  Vnm_tprint( 2, "AVS finite element mesh input not supported!\n");
3674  return VRC_FAILURE;
3675  case VDF_MCSF:
3676  Vnm_tprint(1, "Reading MCSF-format input finite element mesh from %s.\n",
3677  nosh->meshpath[imesh]);
3678  sock = Vio_ctor("FILE", "ASC", VNULL, nosh->meshpath[imesh], "r");
3679  if (sock == VNULL) {
3680  Vnm_print(2, "Problem opening virtual socket %s!\n",
3681  nosh->meshpath[imesh]);
3682  return VRC_FAILURE;
3683  }
3684  if (Vio_accept(sock, 0) < 0) {
3685  Vnm_print(2, "Problem accepting virtual socket %s!\n",
3686  nosh->meshpath[imesh]);
3687  return VRC_FAILURE;
3688  }
3689  break;
3690 
3691  default:
3692  Vnm_tprint( 2, "Invalid data format (%d)!\n",
3693  nosh->meshfmt[imesh]);
3694  return VRC_FAILURE;
3695  }
3696  } else { /* if (feparm->useMesh) */
3697  meshType = VML_DIRICUBE;
3698  for (i=0; i<3; i++) {
3699  center[i] = alist[theMol]->center[i];
3700  length[i] = feparm->glen[i];
3701  }
3702  }
3703 
3704  /* Load the mesh with a particular center and vertex length using the provided input socket */
3705  vrc = Vfetk_loadMesh(fetk[icalc], center, length, meshType, sock);
3706  if (vrc == VRC_FAILURE) {
3707  Vnm_print(2, "Error constructing finite element mesh!\n");
3708  return VRC_FAILURE;
3709  }
3710  //Vnm_redirect(0); // Redirect output to io.mc
3711  Gem_shapeChk(fetk[icalc]->gm); // Traverse simplices and check shapes using the geometry manager.
3712  //Vnm_redirect(1);
3713 
3714  /* Uniformly refine the mesh a bit */
3715  for (j=0; j<2; j++) {
3716  /* AM_* calls below are from the MC package, mc/src/nam/am.c. Note that these calls actually are
3717  * wrappers around Aprx_* functions found in MC as well. */
3718  /* Mark the mesh for needed refinements */
3719  AM_markRefine(fetk[icalc]->am, 0, -1, 0, 0.0);
3720  /* Do actual mesh refinement */
3721  AM_refine(fetk[icalc]->am, 2, 0, feparm->pkey);
3722  //Vnm_redirect(0); // Redirect output to io.mc
3723  Gem_shapeChk(fetk[icalc]->gm); // Traverse simplices and check shapes using the geometry manager.
3724  //Vnm_redirect(1);
3725  }
3726 
3727  /* Setup time statistics */
3728  Vnm_tstop(27, "Setup timer");
3729 
3730  /* Memory statistics */
3731  bytesTotal = Vmem_bytesTotal();
3732  highWater = Vmem_highWaterTotal();
3733 
3734 #ifndef VAPBSQUIET
3735  Vnm_tprint( 1, " Current memory usage: %4.3f MB total, \
3736 %4.3f MB high water\n", (double)(bytesTotal)/(1024.*1024.),
3737  (double)(highWater)/(1024.*1024.));
3738 #endif
3739 
3740  return VRC_SUCCESS;
3741 }
3742 
3743 VPUBLIC void printFEPARM(int icalc,
3744  NOsh *nosh,
3745  FEMparm *feparm,
3746  Vfetk *fetk[NOSH_MAXCALC]
3747  ) {
3748 
3749  Vnm_tprint(1, " Domain size: %g A x %g A x %g A\n",
3750  feparm->glen[0], feparm->glen[1],
3751  feparm->glen[2]);
3752  switch(feparm->ekey) {
3753  case FET_SIMP:
3754  Vnm_tprint(1, " Per-simplex error tolerance: %g\n", feparm->etol);
3755  break;
3756  case FET_GLOB:
3757  Vnm_tprint(1, " Global error tolerance: %g\n", feparm->etol);
3758  break;
3759  case FET_FRAC:
3760  Vnm_tprint(1, " Fraction of simps to refine: %g\n", feparm->etol);
3761  break;
3762  default:
3763  Vnm_tprint(2, "Invalid ekey (%d)!\n", feparm->ekey);
3764  VASSERT(0);
3765  break;
3766  }
3767  switch(feparm->akeyPRE) {
3768  case FRT_UNIF:
3769  Vnm_tprint(1, " Uniform pre-solve refinement.\n");
3770  break;
3771  case FRT_GEOM:
3772  Vnm_tprint(1, " Geometry-based pre-solve refinement.\n");
3773  break;
3774  default:
3775  Vnm_tprint(2, "Invalid akeyPRE (%d)!\n", feparm->akeyPRE);
3776  VASSERT(0);
3777  break;
3778  }
3779  switch(feparm->akeySOLVE) {
3780  case FRT_RESI:
3781  Vnm_tprint(1, " Residual-based a posteriori refinement.\n");
3782  break;
3783  case FRT_DUAL:
3784  Vnm_tprint(1, " Dual-based a posteriori refinement.\n");
3785  break;
3786  case FRT_LOCA:
3787  Vnm_tprint(1, " Local-based a posteriori refinement.\n");
3788  break;
3789  default:
3790  Vnm_tprint(2, "Invalid akeySOLVE (%d)!\n", feparm->akeySOLVE);
3791  break;
3792  }
3793  Vnm_tprint(1, " Refinement of initial mesh to ~%d vertices\n",
3794  feparm->targetNum);
3795  Vnm_tprint(1, " Geometry-based refinment lower bound: %g A\n",
3796  feparm->targetRes);
3797  Vnm_tprint(1, " Maximum number of solve-estimate-refine cycles: %d\n",
3798  feparm->maxsolve);
3799  Vnm_tprint(1, " Maximum number of vertices in mesh: %d\n",
3800  feparm->maxvert);
3801 
3802  /* FOLLOWING IS SOLVER-RELATED; BAIL IF NOT SOLVING */
3803  if (nosh->bogus) return;
3804 #ifdef USE_HB
3805  Vnm_tprint(1, " HB linear solver: AM_hPcg\n");
3806 #else
3807  Vnm_tprint(1, " Non-HB linear solver: ");
3808  switch (fetk[icalc]->lkey) {
3809  case VLT_SLU:
3810  Vnm_print(1, "SLU direct\n");
3811  break;
3812  case VLT_MG:
3813  Vnm_print(1, "multigrid\n");
3814  break;
3815  case VLT_CG:
3816  Vnm_print(1, "conjugate gradient\n");
3817  break;
3818  case VLT_BCG:
3819  Vnm_print(1, "BiCGStab\n");
3820  break;
3821  default:
3822  Vnm_print(1, "???\n");
3823  break;
3824  }
3825 #endif
3826 
3827  Vnm_tprint(1, " Linear solver tol.: %g\n", fetk[icalc]->ltol);
3828  Vnm_tprint(1, " Linear solver max. iters.: %d\n", fetk[icalc]->lmax);
3829  Vnm_tprint(1, " Linear solver preconditioner: ");
3830  switch (fetk[icalc]->lprec) {
3831  case VPT_IDEN:
3832  Vnm_print(1, "identity\n");
3833  break;
3834  case VPT_DIAG:
3835  Vnm_print(1, "diagonal\n");
3836  break;
3837  case VPT_MG:
3838  Vnm_print(1, "multigrid\n");
3839  break;
3840  default:
3841  Vnm_print(1, "???\n");
3842  break;
3843  }
3844  Vnm_tprint(1, " Nonlinear solver: ");
3845  switch (fetk[icalc]->nkey) {
3846  case VNT_NEW:
3847  Vnm_print(1, "newton\n");
3848  break;
3849  case VNT_INC:
3850  Vnm_print(1, "incremental\n");
3851  break;
3852  case VNT_ARC:
3853  Vnm_print(1, "pseudo-arclength\n");
3854  break;
3855  default:
3856  Vnm_print(1, "??? ");
3857  break;
3858  }
3859  Vnm_tprint(1, " Nonlinear solver tol.: %g\n", fetk[icalc]->ntol);
3860  Vnm_tprint(1, " Nonlinear solver max. iters.: %d\n", fetk[icalc]->nmax);
3861  Vnm_tprint(1, " Initial guess: ");
3862  switch (fetk[icalc]->gues) {
3863  case VGT_ZERO:
3864  Vnm_tprint(1, "zero\n");
3865  break;
3866  case VGT_DIRI:
3867  Vnm_tprint(1, "boundary function\n");
3868  break;
3869  case VGT_PREV:
3870  Vnm_tprint(1, "interpolated previous solution\n");
3871  break;
3872  default:
3873  Vnm_tprint(1, "???\n");
3874  break;
3875  }
3876 
3877 }
3878 
3879 VPUBLIC int partFE(int icalc, NOsh *nosh, FEMparm *feparm,
3880  Vfetk *fetk[NOSH_MAXCALC]) {
3881 
3882  Vfetk_setAtomColors(fetk[icalc]);
3883  return 1;
3884 }
3885 
3886 VPUBLIC int preRefineFE(int icalc, /* Calculation index */
3887  FEMparm *feparm, /* FE-specific parameters */
3888  Vfetk *fetk[NOSH_MAXCALC] /* Array of FE solver objects */
3889  ) {
3890 
3891  int nverts,
3892  marked;
3895  /* Based on the refinement type, alert the user to what we're tryng to refine with. */
3896  switch(feparm->akeyPRE) {
3897  case FRT_UNIF:
3898  Vnm_tprint(1, " Commencing uniform refinement to %d verts.\n",
3899  feparm->targetNum);
3900  break;
3901  case FRT_GEOM:
3902  Vnm_tprint(1, " Commencing geometry-based refinement to %d \
3903 verts or %g A resolution.\n", feparm->targetNum, feparm->targetRes);
3904  break;
3905  case FRT_DUAL:
3906  Vnm_tprint(2, "What? You can't do a posteriori error estimation \
3907 before you solve!\n");
3908  VASSERT(0);
3909  break;
3910  case FRT_RESI:
3911  case FRT_LOCA:
3912  default:
3913  VASSERT(0);
3914  break;
3915  }
3916 
3922  Vnm_tprint(1, " Initial mesh has %d vertices\n",
3923  Gem_numVV(fetk[icalc]->gm));
3924 
3925  /* As long as we have simplices marked that need to be refined, run MC's
3926  * AM_markRefine against our data until we hit the error or size tolerance
3927  * for the refinement. */
3928  while (1) {
3929  nverts = Gem_numVV(fetk[icalc]->gm);
3930  if (nverts > feparm->targetNum) {
3931  Vnm_tprint(1, " Hit vertex number limit.\n");
3932  break;
3933  }
3934  marked = AM_markRefine(fetk[icalc]->am, feparm->akeyPRE, -1,
3935  feparm->ekey, feparm->etol);
3936  if (marked == 0) {
3937  Vnm_tprint(1, " Marked 0 simps; hit error/size tolerance.\n");
3938  break;
3939  }
3940  Vnm_tprint(1, " Have %d verts, marked %d. Refining...\n", nverts,
3941  marked);
3942  AM_refine(fetk[icalc]->am, 0, 0, feparm->pkey);
3943  }
3944 
3945  nverts = Gem_numVV(fetk[icalc]->gm);
3946  Vnm_tprint(1, " Done refining; have %d verts.\n", nverts);
3947 
3948  return 1;
3949 }
3950 
3951 
3956 VPUBLIC int solveFE(int icalc,
3957  PBEparm *pbeparm,
3958  FEMparm *feparm,
3959  Vfetk *fetk[NOSH_MAXCALC]
3960  ) {
3961 
3962  int lkeyHB = 3,
3963  meth = 2,
3964  prob = 0,
3965  prec = 0;
3967  if ((pbeparm->pbetype==PBE_NPBE) ||
3968  (pbeparm->pbetype == PBE_NRPBE) ||
3969  (pbeparm->pbetype == PBE_SMPBE)) {
3970 
3971  /* Call MC's nonlinear solver - mc/src/nam/nsolv.c */
3972  AM_nSolve(
3973  fetk[icalc]->am,
3974  fetk[icalc]->nkey,
3975  fetk[icalc]->nmax,
3976  fetk[icalc]->ntol,
3977  meth,
3978  fetk[icalc]->lmax,
3979  fetk[icalc]->ltol,
3980  prec,
3981  fetk[icalc]->gues,
3982  fetk[icalc]->pjac
3983  );
3984  } else if ((pbeparm->pbetype==PBE_LPBE) ||
3985  (pbeparm->pbetype==PBE_LRPBE)) {
3986  /* Note: USEHB is a compile time defined macro. The program flow
3987  is to always take the route using AM_hlSolve when the solver
3988  is linear. D. Gohara 6/29/06
3989  */
3990 #ifdef USE_HB
3991  Vnm_print(2, "SORRY! DON'T USE HB!!!\n");
3992  VASSERT(0);
3993 
3994  /* Call MC's hierarchical linear solver - mc/src/nam/lsolv.c */
3995  AM_hlSolve(fetk[icalc]->am, meth, lkeyHB, fetk[icalc]->lmax,
3996  fetk[icalc]->ltol, fetk[icalc]->gues, fetk[icalc]->pjac);
3997 #else
3998 
3999  /* Call MC's linear solver - mc/src/nam/lsolv.c */
4000  AM_lSolve(
4001  fetk[icalc]->am,
4002  prob,
4003  meth,
4004  fetk[icalc]->lmax,
4005  fetk[icalc]->ltol,
4006  prec,
4007  fetk[icalc]->gues,
4008  fetk[icalc]->pjac
4009  );
4010 #endif
4011  }
4012 
4013  return 1;
4014 }
4015 
4019 VPUBLIC int energyFE(NOsh *nosh,
4020  int icalc,
4021  Vfetk *fetk[NOSH_MAXCALC],
4022  int *nenergy,
4023  double *totEnergy,
4024  double *qfEnergy,
4025  double *qmEnergy,
4026  double *dielEnergy
4027  ) {
4028 
4029  FEMparm *feparm = nosh->calc[icalc]->femparm;
4030  PBEparm *pbeparm = nosh->calc[icalc]->pbeparm;
4032  *nenergy = 1;
4033  *totEnergy = 0;
4034 
4041  if (nosh->bogus == 0) {
4042  if ((pbeparm->pbetype==PBE_NPBE) ||
4043  (pbeparm->pbetype==PBE_NRPBE) ||
4044  (pbeparm->pbetype == PBE_SMPBE)) {
4045  *totEnergy = Vfetk_energy(fetk[icalc], -1, 1); /* Last parameter indicates NPBE */
4046  } else if ((pbeparm->pbetype==PBE_LPBE) ||
4047  (pbeparm->pbetype==PBE_LRPBE)) {
4048  *totEnergy = Vfetk_energy(fetk[icalc], -1, 0); /* Last parameter indicates LPBE */
4049  } else {
4050  VASSERT(0);
4051  }
4052 
4053 #ifndef VAPBSQUIET
4054  Vnm_tprint(1, " Total electrostatic energy = %1.12E kJ/mol\n",
4055  Vunit_kb*pbeparm->temp*(1e-3)*Vunit_Na*(*totEnergy));
4056  fflush(stdout);
4057 #endif
4058  }
4059 
4060  if (pbeparm->calcenergy == PCE_COMPS) {
4061  Vnm_tprint(2, "Error! Verbose energy evaluation not available for FEM yet!\n");
4062  Vnm_tprint(2, "E-mail nathan.baker@pnl.gov if you want this.\n");
4063  *qfEnergy = 0;
4064  *qmEnergy = 0;
4065  *dielEnergy = 0;
4066  } else {
4067  *nenergy = 0;
4068  }
4069  return 1;
4070 }
4071 
4079 VPUBLIC int postRefineFE(int icalc,
4080  FEMparm *feparm,
4081  Vfetk *fetk[NOSH_MAXCALC]
4082  ) {
4083 
4084  int nverts,
4085  marked;
4087  nverts = Gem_numVV(fetk[icalc]->gm);
4088  if (nverts > feparm->maxvert) {
4089  Vnm_tprint(1, " Current number of vertices (%d) exceeds max (%d)!\n",
4090  nverts, feparm->maxvert);
4091  return 0;
4092  }
4093  Vnm_tprint(1, " Mesh currently has %d vertices\n", nverts);
4094 
4095  switch(feparm->akeySOLVE) {
4096  case FRT_UNIF:
4097  Vnm_tprint(1, " Commencing uniform refinement.\n");
4098  break;
4099  case FRT_GEOM:
4100  Vnm_tprint(1, " Commencing geometry-based refinement.\n");
4101  break;
4102  case FRT_RESI:
4103  Vnm_tprint(1, " Commencing residual-based refinement.\n");
4104  break;
4105  case FRT_DUAL:
4106  Vnm_tprint(1, " Commencing dual problem-based refinement.\n");
4107  break;
4108  case FRT_LOCA:
4109  Vnm_tprint(1, " Commencing local-based refinement.\n.");
4110  break;
4111  default:
4112  Vnm_tprint(2, " Error -- unknown refinement type (%d)!\n",
4113  feparm->akeySOLVE);
4114  return 0;
4115  break;
4116  }
4117 
4118  /* Run MC's refinement algorithm */
4119  marked = AM_markRefine(fetk[icalc]->am, feparm->akeySOLVE, -1,
4120  feparm->ekey, feparm->etol);
4121  if (marked == 0) {
4122  Vnm_tprint(1, " Marked 0 simps; hit error/size tolerance.\n");
4123  return 0;
4124  }
4125  Vnm_tprint(1, " Have %d verts, marked %d. Refining...\n", nverts,
4126  marked);
4127  AM_refine(fetk[icalc]->am, 0, 0, feparm->pkey);
4128  nverts = Gem_numVV(fetk[icalc]->gm);
4129  Vnm_tprint(1, " Done refining; have %d verts.\n", nverts);
4130  //Vnm_redirect(0); // Redirect output to io.mc
4131  Gem_shapeChk(fetk[icalc]->gm); // Traverse simplices and check shapes using the geometry manager.
4132  //Vnm_redirect(1);
4133 
4134  return 1;
4135 }
4136 
4140 VPUBLIC int writedataFE(int rank,
4141  NOsh *nosh,
4142  PBEparm *pbeparm,
4143  Vfetk *fetk
4144  ) {
4145 
4146  char writestem[VMAX_ARGLEN];
4147  char outpath[VMAX_ARGLEN];
4148  int i,
4149  writeit;
4150  AM *am;
4151  Bvec *vec;
4153  if (nosh->bogus) return 1;
4154 
4155  am = fetk->am;
4156  vec = am->w0;
4157 
4158  for (i=0; i<pbeparm->numwrite; i++) {
4159 
4160  writeit = 1;
4161 
4162  switch (pbeparm->writetype[i]) {
4163 
4164  case VDT_CHARGE:
4165 
4166  Vnm_tprint(2, " Sorry; can't write charge distribution for FEM!\n");
4167  writeit = 0;
4168  break;
4169 
4170  case VDT_POT:
4171 
4172  Vnm_tprint(1, " Writing potential to ");
4173  Vfetk_fillArray(fetk, vec, VDT_POT);
4174  break;
4175 
4176  case VDT_SMOL:
4177 
4178  Vnm_tprint(1, " Writing molecular accessibility to ");
4179  Vfetk_fillArray(fetk, vec, VDT_SMOL);
4180  break;
4181 
4182  case VDT_SSPL:
4183 
4184  Vnm_tprint(1, " Writing spline-based accessibility to ");
4185  Vfetk_fillArray(fetk, vec, VDT_SSPL);
4186  break;
4187 
4188  case VDT_VDW:
4189 
4190  Vnm_tprint(1, " Writing van der Waals accessibility to ");
4191  Vfetk_fillArray(fetk, vec, VDT_VDW);
4192  break;
4193 
4194  case VDT_IVDW:
4195 
4196  Vnm_tprint(1, " Writing ion accessibility to ");
4197  Vfetk_fillArray(fetk, vec, VDT_IVDW);
4198  break;
4199 
4200  case VDT_LAP:
4201 
4202  Vnm_tprint(2, " Sorry; can't write charge distribution for FEM!\n");
4203  writeit = 0;
4204  break;
4205 
4206  case VDT_EDENS:
4207 
4208  Vnm_tprint(2, " Sorry; can't write energy density for FEM!\n");
4209  writeit = 0;
4210  break;
4211 
4212  case VDT_NDENS:
4213 
4214  Vnm_tprint(1, " Writing number density to ");
4215  Vfetk_fillArray(fetk, vec, VDT_NDENS);
4216  break;
4217 
4218  case VDT_QDENS:
4219 
4220  Vnm_tprint(1, " Writing charge density to ");
4221  Vfetk_fillArray(fetk, vec, VDT_QDENS);
4222  break;
4223 
4224  case VDT_DIELX:
4225 
4226  Vnm_tprint(2, " Sorry; can't write x-shifted dielectric map for FEM!\n");
4227  writeit = 0;
4228  break;
4229 
4230  case VDT_DIELY:
4231 
4232  Vnm_tprint(2, " Sorry; can't write y-shifted dielectric map for FEM!\n");
4233  writeit = 0;
4234  break;
4235 
4236  case VDT_DIELZ:
4237 
4238  Vnm_tprint(2, " Sorry; can't write z-shifted dielectric map for FEM!\n");
4239  writeit = 0;
4240  break;
4241 
4242  case VDT_KAPPA:
4243 
4244  Vnm_tprint(1, " Sorry; can't write kappa map for FEM!\n");
4245  writeit = 0;
4246  break;
4247 
4248  case VDT_ATOMPOT:
4249 
4250  Vnm_tprint(1, " Sorry; can't write atom potentials for FEM!\n");
4251  writeit = 0;
4252  break;
4253 
4254  default:
4255 
4256  Vnm_tprint(2, "Invalid data type for writing!\n");
4257  writeit = 0;
4258  return 0;
4259  }
4260 
4261  if (!writeit) return 0;
4262 
4263 
4264 #ifdef HAVE_MPI_H
4265  sprintf(writestem, "%s-PE%d", pbeparm->writestem[i], rank);
4266 #else
4267  if(nosh->ispara){
4268  sprintf(writestem, "%s-PE%d", pbeparm->writestem[i],nosh->proc_rank);
4269  }else{
4270  sprintf(writestem, "%s", pbeparm->writestem[i]);
4271  }
4272 #endif
4273 
4274  switch (pbeparm->writefmt[i]) {
4275 
4276  case VDF_DX:
4277  sprintf(outpath, "%s.%s", writestem, "dx");
4278  Vnm_tprint(1, "%s\n", outpath);
4279  Vfetk_write(fetk, "FILE", "ASC", VNULL, outpath, vec, VDF_DX);
4280  break;
4281 
4282  case VDF_AVS:
4283  sprintf(outpath, "%s.%s", writestem, "ucd");
4284  Vnm_tprint(1, "%s\n", outpath);
4285  Vfetk_write(fetk, "FILE", "ASC", VNULL, outpath, vec, VDF_AVS);
4286  break;
4287 
4288  case VDF_UHBD:
4289  Vnm_tprint(2, "UHBD format not supported for FEM!\n");
4290  break;
4291 
4292  case VDF_MCSF:
4293  Vnm_tprint(2, "MCSF format not supported yet!\n");
4294  break;
4295 
4296  default:
4297  Vnm_tprint(2, "Bogus data format (%d)!\n",
4298  pbeparm->writefmt[i]);
4299  break;
4300  }
4301 
4302  }
4303 
4304  return 1;
4305 }
4306 #endif /* ifdef HAVE_MCX_H */
4307 
4308 VPUBLIC int initAPOL(NOsh *nosh,
4309  Vmem *mem,
4310  Vparam *param,
4311  APOLparm *apolparm,
4312  int *nforce,
4313  AtomForce **atomForce,
4314  Valist *alist
4315  ) {
4316  int i,
4317  natoms,
4318  len,
4319  inhash[3],
4320  rc = 0;
4322  time_t ts;
4323  Vclist *clist = VNULL;
4324  Vacc *acc = VNULL;
4325  Vatom *atom = VNULL;
4326  Vparam_AtomData *atomData = VNULL;
4328  double sasa,
4329  sav,
4330  nhash[3],
4331  sradPad,
4332  x,
4333  y,
4334  z,
4335  atomRadius,
4336  srad,
4337  *atomsasa,
4338  *atomwcaEnergy,
4339  energy = 0.0,
4340  dist,
4341  charge,
4342  xmin,
4343  xmax,
4344  ymin,
4345  ymax,
4346  zmin,
4347  zmax,
4348  disp[3],
4349  center[3],
4350  soluteXlen,
4351  soluteYlen,
4352  soluteZlen;
4354  atomsasa = (double *)Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
4355  atomwcaEnergy = (double *)Vmem_malloc(VNULL, Valist_getNumberAtoms(alist), sizeof(double));
4356 
4357  /* Determine solute length and charge*/
4358  atom = Valist_getAtom(alist, 0);
4359  xmin = Vatom_getPosition(atom)[0];
4360  xmax = Vatom_getPosition(atom)[0];
4361  ymin = Vatom_getPosition(atom)[1];
4362  ymax = Vatom_getPosition(atom)[1];
4363  zmin = Vatom_getPosition(atom)[2];
4364  zmax = Vatom_getPosition(atom)[2];
4365  charge = 0;
4366  natoms = Valist_getNumberAtoms(alist);
4367 
4368  for (i=0; i < natoms; i++) {
4369  atom = Valist_getAtom(alist, i);
4370  atomRadius = Vatom_getRadius(atom);
4371  x = Vatom_getPosition(atom)[0];
4372  y = Vatom_getPosition(atom)[1];
4373  z = Vatom_getPosition(atom)[2];
4374  if ((x+atomRadius) > xmax) xmax = x + atomRadius;
4375  if ((x-atomRadius) < xmin) xmin = x - atomRadius;
4376  if ((y+atomRadius) > ymax) ymax = y + atomRadius;
4377  if ((y-atomRadius) < ymin) ymin = y - atomRadius;
4378  if ((z+atomRadius) > zmax) zmax = z + atomRadius;
4379  if ((z-atomRadius) < zmin) zmin = z - atomRadius;
4380  disp[0] = (x - center[0]);
4381  disp[1] = (y - center[1]);
4382  disp[2] = (z - center[2]);
4383  dist = (disp[0]*disp[0]) + (disp[1]*disp[1]) + (disp[2]*disp[2]);
4384  dist = VSQRT(dist) + atomRadius;
4385  charge += Vatom_getCharge(Valist_getAtom(alist, i));
4386  }
4387  soluteXlen = xmax - xmin;
4388  soluteYlen = ymax - ymin;
4389  soluteZlen = zmax - zmin;
4390 
4391  /* Set up the hash table for the cell list */
4392  Vnm_print(0, "APOL: Setting up hash table and accessibility object...\n");
4393  nhash[0] = soluteXlen/0.5;
4394  nhash[1] = soluteYlen/0.5;
4395  nhash[2] = soluteZlen/0.5;
4396  for (i=0; i<3; i++) inhash[i] = (int)(nhash[i]);
4397 
4398  for (i=0;i<3;i++){
4399  if (inhash[i] < 3) inhash[i] = 3;
4400  if (inhash[i] > MAX_HASH_DIM) inhash[i] = MAX_HASH_DIM;
4401  }
4402 
4403  /* Pad the radius by 2x the maximum displacement value */
4404  srad = apolparm->srad;
4405  sradPad = srad + (2*apolparm->dpos);
4406  clist = Vclist_ctor(alist, sradPad , inhash, CLIST_AUTO_DOMAIN,
4407  VNULL, VNULL);
4408  acc = Vacc_ctor(alist, clist, apolparm->sdens);
4409 
4410  /* Get WAT (water) LJ parameters from Vparam object */
4411  if (param == VNULL && (apolparm->bconc != 0.0)) {
4412  Vnm_tprint(2, "initAPOL: Got NULL Vparam object!\n");
4413  Vnm_tprint(2, "initAPOL: You are performing an apolar calculation with the van der Waals integral term,\n");
4414  Vnm_tprint(2, "initAPOL: this term requires van der Waals parameters which are not available from the \n");
4415  Vnm_tprint(2, "initAPOL: PQR file. Therefore, you need to supply a parameter file with the parm keyword,\n");
4416  Vnm_tprint(2, "initAPOL: for example,\n");
4417  Vnm_tprint(2, "initAPOL: read parm flat amber94.dat end\n");
4418  Vnm_tprint(2, "initAPOL: where the relevant parameter files can be found in apbs/tools/conversion/param/vparam.\n");
4419  return VRC_FAILURE;
4420  }
4421 
4422  if (apolparm->bconc != 0.0){
4423  atomData = Vparam_getAtomData(param, "WAT", "OW");
4424  if (atomData == VNULL) atomData = Vparam_getAtomData(param, "WAT", "O");
4425  if (atomData == VNULL) {
4426  Vnm_tprint(2, "initAPOL: Couldn't find parameters for WAT OW or WAT O!\n");
4427  Vnm_tprint(2, "initAPOL: These parameters must be present in your file\n");
4428  Vnm_tprint(2, "initAPOL: for apolar calculations.\n");
4429  return VRC_FAILURE;
4430  }
4431  apolparm->watepsilon = atomData->epsilon;
4432  apolparm->watsigma = atomData->radius;
4433  apolparm->setwat = 1;
4434  }
4435 
4436  /* Calculate Energy and Forces */
4437  if(apolparm->calcforce) {
4438  rc = forceAPOL(acc, mem, apolparm, nforce, atomForce, alist, clist);
4439  if(rc == VRC_FAILURE) {
4440  Vnm_print(2, "Error in apolar force calculation!\n");
4441  return VRC_FAILURE;
4442  }
4443  }
4444 
4445  /* Get the SAV and SAS */
4446  sasa = 0.0;
4447  sav = 0.0;
4448 
4449  if (apolparm->calcenergy) {
4450  len = Valist_getNumberAtoms(alist);
4451 
4452  if (VABS(apolparm->gamma) > VSMALL) {
4453  /* Total Solvent Accessible Surface Area (SASA) */
4454  apolparm->sasa = Vacc_totalSASA(acc, srad);
4455  /* SASA for each atom */
4456  for (i = 0; i < len; i++) {
4457  atom = Valist_getAtom(alist, i);
4458  atomsasa[i] = Vacc_atomSASA(acc, srad, atom);
4459  }
4460  } else {
4461  /* Total Solvent Accessible Surface Area (SASA) set to zero */
4462  apolparm->sasa = 0.0;
4463  /* SASA for each atom set to zero*/
4464  for (i = 0; i < len; i++) {
4465  atomsasa[i] = 0.0;
4466  }
4467  }
4468 
4469  /* Inflated van der Waals accessibility */
4470  if (VABS(apolparm->press) > VSMALL){
4471  apolparm->sav = Vacc_totalSAV(acc, clist, apolparm, srad);
4472  } else {
4473  apolparm->sav = 0.0;
4474  }
4475 
4476  /* wcaEnergy integral code */
4477  if (VABS(apolparm->bconc) > VSMALL) {
4478  /* wcaEnergy for each atom */
4479  for (i = 0; i < len; i++) {
4480  rc = Vacc_wcaEnergyAtom(acc, apolparm, alist, clist, i, &energy);
4481  if (rc == 0) {
4482  Vnm_print(2, "Error in apolar energy calculation!\n");
4483  return 0;
4484  }
4485  atomwcaEnergy[i] = energy;
4486  }
4487  /* Total WCA Energy */
4488  rc = Vacc_wcaEnergy(acc, apolparm, alist, clist);
4489  if (rc == 0) {
4490  Vnm_print(2, "Error in apolar energy calculation!\n");
4491  return 0;
4492  }
4493  } else {
4494  apolparm->wcaEnergy = 0.0;
4495  }
4496  energyAPOL(apolparm, apolparm->sasa, apolparm->sav, atomsasa, atomwcaEnergy, Valist_getNumberAtoms(alist));
4497  }
4498 
4499  Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double), (void **)&(atomsasa));
4500  Vmem_free(VNULL, Valist_getNumberAtoms(alist), sizeof(double), (void **)&(atomwcaEnergy));
4501  Vclist_dtor(&clist);
4502  Vacc_dtor(&acc);
4503 
4504  return VRC_SUCCESS;
4505 }
4506 
4507 VPUBLIC int energyAPOL(APOLparm *apolparm,
4508  double sasa,
4509  double sav,
4510  double atomsasa[],
4511  double atomwcaEnergy[],
4512  int numatoms
4513  ){
4514 
4515  double energy = 0.0;
4516  int i = 0;
4517 
4518 #ifndef VAPBSQUIET
4519  Vnm_print(1,"\nSolvent Accessible Surface Area (SASA) for each atom:\n");
4520  for (i = 0; i < numatoms; i++) {
4521  Vnm_print(1," SASA for atom %i: %1.12E\n", i, atomsasa[i]);
4522  }
4523 
4524  Vnm_print(1,"\nTotal solvent accessible surface area: %g A^2\n",sasa);
4525 #endif
4526 
4527  switch(apolparm->calcenergy){
4528  case ACE_NO:
4529  break;
4530  case ACE_COMPS:
4531  Vnm_print(1,"energyAPOL: Cannot calculate component energy, skipping.\n");
4532  break;
4533  case ACE_TOTAL:
4534  energy = (apolparm->gamma*sasa) + (apolparm->press*sav)
4535  + (apolparm->wcaEnergy);
4536 #ifndef VAPBSQUIET
4537  Vnm_print(1,"\nSurface tension*area energies (gamma * SASA) for each atom:\n");
4538  for (i = 0; i < numatoms; i++) {
4539  Vnm_print(1," Surface tension*area energy for atom %i: %1.12E\n", i, apolparm->gamma*atomsasa[i]);
4540  }
4541 
4542  Vnm_print(1,"\nTotal surface tension energy: %g kJ/mol\n", apolparm->gamma*sasa);
4543  Vnm_print(1,"\nTotal solvent accessible volume: %g A^3\n", sav);
4544  Vnm_print(1,"\nTotal pressure*volume energy: %g kJ/mol\n", apolparm->press*sav);
4545  Vnm_print(1,"\nWCA dispersion Energies for each atom:\n");
4546  for (i = 0; i < numatoms; i++) {
4547  Vnm_print(1," WCA energy for atom %i: %1.12E\n", i, atomwcaEnergy[i]);
4548  }
4549 
4550  Vnm_print(1,"\nTotal WCA energy: %g kJ/mol\n", (apolparm->wcaEnergy));
4551  Vnm_print(1,"\nTotal non-polar energy = %1.12E kJ/mol\n", energy);
4552 #endif
4553  break;
4554  default:
4555  Vnm_print(2,"energyAPOL: Error in energyAPOL. Unknown option.\n");
4556  break;
4557  }
4558 
4559  return VRC_SUCCESS;
4560 }
4561 
4562 VPUBLIC int forceAPOL(Vacc *acc,
4563  Vmem *mem,
4564  APOLparm *apolparm,
4565  int *nforce,
4566  AtomForce **atomForce,
4567  Valist *alist,
4568  Vclist *clist
4569  ) {
4570  time_t ts, ts_main, ts_sub;
4571  int i,
4572  j,
4573  natom;
4574 
4575  double srad, /* Probe radius */
4576  xF,
4577  yF,
4578  zF, /* Individual forces */
4579  press,
4580  gamma,
4581  offset,
4582  bconc,
4583  dSASA[3],
4584  dSAV[3],
4585  force[3],
4586  *apos;
4587 
4588  Vatom *atom = VNULL;
4589  ts_main = clock();
4590 
4591  srad = apolparm->srad;
4592  press = apolparm->press;
4593  gamma = apolparm->gamma;
4594  offset = apolparm->dpos;
4595  bconc = apolparm->bconc;
4596 
4597  natom = Valist_getNumberAtoms(alist);
4598 
4599  /* Check to see if we need to build the surface */
4600  Vnm_print(0, "forceAPOL: Trying atom surf...\n");
4601  ts = clock();
4602  if (acc->surf == VNULL) {
4603  acc->surf = (VaccSurf**)Vmem_malloc(acc->mem, natom, sizeof(VaccSurf *));
4604  for (i=0; i<natom; i++) {
4605  atom = Valist_getAtom(acc->alist, i);
4606  //apos = Vatom_getPosition(atom); // apos never referenced? - Peter
4607  /* NOTE: RIGHT NOW WE DO THIS FOR THE ENTIRE MOLECULE WHICH IS
4608  * INCREDIBLY INEFFICIENT, PARTICULARLY DURING FOCUSING!!! */
4609  acc->surf[i] = Vacc_atomSurf(acc, atom, acc->refSphere, srad);
4610  }
4611  }
4612  Vnm_print(0, "forceAPOL: atom surf: Time elapsed: %f\n", ((double)clock() - ts) / CLOCKS_PER_SEC);
4613 
4614  if(apolparm->calcforce == ACF_TOTAL){
4615  Vnm_print(0, "forceAPOL: calcforce == ACF_TOTAL\n");
4616  ts = clock();
4617 
4618  *nforce = 1;
4619  if(*atomForce != VNULL){
4620  Vmem_free(mem,*nforce,sizeof(AtomForce), (void **)atomForce);
4621  }
4622 
4623  *atomForce = (AtomForce *)Vmem_malloc(mem, *nforce,
4624  sizeof(AtomForce));
4625 
4626  /* Clear out force arrays */
4627  for (j=0; j<3; j++) {
4628  (*atomForce)[0].sasaForce[j] = 0.0;
4629  (*atomForce)[0].savForce[j] = 0.0;
4630  (*atomForce)[0].wcaForce[j] = 0.0;
4631  }
4632 
4633  // problem block
4634  for (i=0; i<natom; i++) {
4635  atom = Valist_getAtom(alist, i);
4636 
4637  for(j=0;j<3;j++){
4638  dSASA[j] = 0.0;
4639  dSAV[j] = 0.0;
4640  force[j] = 0.0;
4641  }
4642 
4643  if(VABS(gamma) > VSMALL) {
4644  Vacc_atomdSASA(acc, offset, srad, atom, dSASA);
4645  }
4646  if(VABS(press) > VSMALL) {
4647  Vacc_atomdSAV(acc, srad, atom, dSAV);
4648  }
4649  if(VABS(bconc) > VSMALL) {
4650  Vacc_wcaForceAtom(acc, apolparm, clist, atom, force);
4651  }
4652 
4653  for(j=0;j<3;j++){
4654  (*atomForce)[0].sasaForce[j] += dSASA[j];
4655  (*atomForce)[0].savForce[j] += dSAV[j];
4656  (*atomForce)[0].wcaForce[j] += force[j];
4657  }
4658  }
4659  // end block
4660 
4661  Vnm_tprint( 1, " Printing net forces (kJ/mol/A)\n");
4662  Vnm_tprint( 1, " Legend:\n");
4663  Vnm_tprint( 1, " sasa -- SASA force\n");
4664  Vnm_tprint( 1, " sav -- SAV force\n");
4665  Vnm_tprint( 1, " wca -- WCA force\n\n");
4666 
4667  Vnm_tprint( 1, " sasa %4.3e %4.3e %4.3e\n",
4668  (*atomForce)[0].sasaForce[0],
4669  (*atomForce)[0].sasaForce[1],
4670  (*atomForce)[0].sasaForce[2]);
4671  Vnm_tprint( 1, " sav %4.3e %4.3e %4.3e\n",
4672  (*atomForce)[0].savForce[0],
4673  (*atomForce)[0].savForce[1],
4674  (*atomForce)[0].savForce[2]);
4675  Vnm_tprint( 1, " wca %4.3e %4.3e %4.3e\n",
4676  (*atomForce)[0].wcaForce[0],
4677  (*atomForce)[0].wcaForce[1],
4678  (*atomForce)[0].wcaForce[2]);
4679 
4680  Vnm_print(0, "forceAPOL: calcforce == ACF_TOTAL: %f\n", ((double)clock() - ts) / CLOCKS_PER_SEC);
4681  } else if (apolparm->calcforce == ACF_COMPS ){
4682  *nforce = Valist_getNumberAtoms(alist);
4683  if(*atomForce == VNULL){
4684  *atomForce = (AtomForce *)Vmem_malloc(mem, *nforce,
4685  sizeof(AtomForce));
4686  }else{
4687  Vmem_free(mem,*nforce,sizeof(AtomForce), (void **)atomForce);
4688  *atomForce = (AtomForce *)Vmem_malloc(mem, *nforce,
4689  sizeof(AtomForce));
4690  }
4691 
4692 #ifndef VAPBSQUIET
4693  Vnm_tprint( 1, " Printing per atom forces (kJ/mol/A)\n");
4694  Vnm_tprint( 1, " Legend:\n");
4695  Vnm_tprint( 1, " tot n -- Total force for atom n\n");
4696  Vnm_tprint( 1, " sasa n -- SASA force for atom n\n");
4697  Vnm_tprint( 1, " sav n -- SAV force for atom n\n");
4698  Vnm_tprint( 1, " wca n -- WCA force for atom n\n\n");
4699 
4700  Vnm_tprint( 1, " gamma %f\n" \
4701  " pressure %f\n" \
4702  " bconc %f \n\n",
4703  gamma,press,bconc);
4704 #endif
4705 
4706  for (i=0; i<natom; i++) {
4707  atom = Valist_getAtom(alist, i);
4708 
4709  for(j=0;j<3;j++){
4710  dSASA[j] = 0.0;
4711  dSAV[j] = 0.0;
4712  force[j] = 0.0;
4713  }
4714 
4715  /* Clear out force arrays */
4716  for (j=0; j<3; j++) {
4717  (*atomForce)[i].sasaForce[j] = 0.0;
4718  (*atomForce)[i].savForce[j] = 0.0;
4719  (*atomForce)[i].wcaForce[j] = 0.0;
4720  }
4721 
4722  if(VABS(gamma) > VSMALL) Vacc_atomdSASA(acc, offset, srad, atom, dSASA);
4723  if(VABS(press) > VSMALL) Vacc_atomdSAV(acc, srad, atom, dSAV);
4724  if(VABS(bconc) > VSMALL) Vacc_wcaForceAtom(acc,apolparm,clist,atom,force);
4725 
4726  xF = -((gamma*dSASA[0]) + (press*dSAV[0]) + (bconc*force[0]));
4727  yF = -((gamma*dSASA[1]) + (press*dSAV[1]) + (bconc*force[1]));
4728  zF = -((gamma*dSASA[2]) + (press*dSAV[2]) + (bconc*force[2]));
4729 
4730  for(j=0;j<3;j++){
4731  (*atomForce)[i].sasaForce[j] += dSASA[j];
4732  (*atomForce)[i].savForce[j] += dSAV[j];
4733  (*atomForce)[i].wcaForce[j] += force[j];
4734  }
4735 
4736 #ifndef VAPBSQUIET
4737  Vnm_print( 1, " tot %i %4.3e %4.3e %4.3e\n",
4738  i,
4739  xF,
4740  yF,
4741  zF);
4742  Vnm_print( 1, " sasa %i %4.3e %4.3e %4.3e\n",
4743  i,
4744  (*atomForce)[i].sasaForce[0],
4745  (*atomForce)[i].sasaForce[1],
4746  (*atomForce)[i].sasaForce[2]);
4747  Vnm_print( 1, " sav %i %4.3e %4.3e %4.3e\n",
4748  i,
4749  (*atomForce)[i].savForce[0],
4750  (*atomForce)[i].savForce[1],
4751  (*atomForce)[i].savForce[2]);
4752  Vnm_print( 1, " wca %i %4.3e %4.3e %4.3e\n",
4753  i,
4754  (*atomForce)[i].wcaForce[0],
4755  (*atomForce)[i].wcaForce[1],
4756  (*atomForce)[i].wcaForce[2]);
4757 #endif
4758 
4759  }
4760  } else *nforce = 0;
4761 
4762 #ifndef VAPBSQUIET
4763  Vnm_print(1,"\n");
4764 #endif
4765 
4766  Vnm_print(0, "forceAPOL: Time elapsed: %f\n", ((double)clock() - ts_main) / CLOCKS_PER_SEC);
4767  return VRC_SUCCESS;
4768 }
4769 
4770 
Definition: vfetk.h:89
Definition: vhal.h:273
double glen[3]
Definition: femparm.h:140
double sav
Definition: apolparm.h:176
VPUBLIC Vacc * Vacc_ctor(Valist *alist, Vclist *clist, double surf_density)
Construct the accessibility object.
Definition: vacc.c:132
double zmax
Definition: vpmgp.h:197
VPUBLIC double Vpmg_qmEnergy(Vpmg *thee, int extFlag)
Get the "mobile charge" contribution to the electrostatic energy.
Definition: vpmg.c:1391
double targetRes
Definition: femparm.h:160
Definition: vhal.h:316
Definition: vfetk.h:122
Definition: vhal.h:217
MGparm * mgparm
Definition: nosh.h:165
VPUBLIC double Vpmg_qfEnergy(Vpmg *thee, int extFlag)
Get the "fixed charge" contribution to the electrostatic energy.
Definition: vpmg.c:1692
Vdata_Format meshfmt[NOSH_MAXMOL]
Definition: nosh.h:246
double bconc
Definition: apolparm.h:139
VPUBLIC void Vgrid_writeDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the data in OpenDX grid format.
Definition: vgrid.c:1001
int meshID
Definition: femparm.h:174
Vmem * vmem
Definition: vpmg.h:118
double hx
Definition: vpmgp.h:87
int apol2calc[NOSH_MAXCALC]
Definition: nosh.h:217
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:1592
VPUBLIC Vparam * loadParameter(NOsh *nosh)
Loads and returns parameter object.
Definition: routines.c:60
double ymax
Definition: vpmgp.h:196
int targetNum
Definition: femparm.h:155
VPUBLIC void Vfetk_setAtomColors(Vfetk *thee)
Transfer color (partition ID) information frmo a partitioned mesh to the atoms.
Definition: vfetk.c:849
int potMapID
Definition: pbeparm.h:129
VPUBLIC void Vgrid_writeUHBD(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out the data in UHBD grid format.
Definition: vgrid.c:1251
VPUBLIC Vfetk * Vfetk_ctor(Vpbe *pbe, Vhal_PBEType type)
Constructor for Vfetk object.
Definition: vfetk.c:532
int gotparm
Definition: nosh.h:224
int useKappaMap
Definition: pbeparm.h:124
double epsilon
Definition: vparam.h:97
VPUBLIC double Vacc_atomSASA(Vacc *thee, double radius, Vatom *atom)
Return the atomic solvent accessible surface area (SASA)
Definition: vacc.c:780
#define NOSH_MAXMOL
Maximum number of molecules in a run.
Definition: nosh.h:79
double temp
Definition: apolparm.h:160
int chargeMapID
Definition: pbeparm.h:133
Contains public data members for Vpbe class/module.
Definition: vpbe.h:84
Oracle for solvent- and ion-accessibility around a biomolecule.
Definition: vacc.h:108
Vdata_Format dielfmt[NOSH_MAXMOL]
Definition: nosh.h:234
VPUBLIC Vparam * Vparam_ctor()
Construct the object.
Definition: vparam.c:181
VPUBLIC int loadChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the charge maps given in NOsh into grid objects.
Definition: routines.c:770
double Lmem
Definition: pbeparm.h:176
double ymin
Definition: vpmgp.h:193
VPUBLIC void Vpmgp_dtor(Vpmgp **thee)
Object destructor.
Definition: vpmgp.c:178
VPUBLIC void Vgrid_dtor(Vgrid **thee)
Object destructor.
Definition: vgrid.c:143
double hzed
Definition: vpmgp.h:89
double sdens
Definition: apolparm.h:142
VPUBLIC void Vacc_dtor(Vacc **thee)
Destroy object.
Definition: vacc.c:245
char potpath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:239
double xmin
Definition: vpmgp.h:192
VaccSurf ** surf
Definition: vacc.h:117
VPUBLIC int energyAPOL(APOLparm *apolparm, double sasa, double sav, double atomsasa[], double atomwcaEnergy[], int numatoms)
Calculate non-polar energies.
Definition: routines.c:4507
NOsh_PrintType printwhat[NOSH_MAXPRINT]
Definition: nosh.h:248
int molid
Definition: pbeparm.h:119
Electrostatic potential oracle for Cartesian mesh data.
Definition: vgrid.h:81
double smvolume
Definition: pbeparm.h:162
int nprint
Definition: nosh.h:247
VPUBLIC int Vgrid_readGZ(Vgrid *thee, const char *fname)
Read in OpenDX data in GZIP format.
Definition: vgrid.c:451
VPUBLIC int writedataFlat(NOsh *nosh, Vcom *com, const char *fname, double totEnergy[NOSH_MAXCALC], double qfEnergy[NOSH_MAXCALC], double qmEnergy[NOSH_MAXCALC], double dielEnergy[NOSH_MAXCALC], int nenergy[NOSH_MAXCALC], double *atomEnergy[NOSH_MAXCALC], int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Write out information to a flat file.
Definition: routines.c:1741
Definition: nosh.h:99
Vdata_Format writefmt[PBEPARM_MAXWRITE]
Definition: pbeparm.h:189
int pkey
Definition: femparm.h:170
VaccSurf * refSphere
Definition: vacc.h:116
char elecname[NOSH_MAXCALC][VMAX_ARGLEN]
Definition: nosh.h:255
double ibForce[3]
Definition: routines.h:84
double hy
Definition: vpmgp.h:88
double swin
Definition: pbeparm.h:154
double xcent
Definition: vpmgp.h:118
VPUBLIC int postRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Estimate error, mark mesh, and refine mesh after solve.
Definition: routines.c:4079
double zcent
Definition: vpmgp.h:120
NOsh_calc * calc[NOSH_MAXCALC]
Definition: nosh.h:185
double srad
Definition: apolparm.h:154
Vhal_PBEType pbetype
Definition: pbeparm.h:134
Definition: nosh.h:98
Contains public data members for Vpmg class/module.
Definition: vpmg.h:116
double * u
Definition: vpmg.h:147
VPUBLIC int Valist_getNumberAtoms(Valist *thee)
Get number of atoms in the list.
Definition: valist.c:105
VPUBLIC int solveMG(NOsh *nosh, Vpmg *pmg, MGparm_CalcType type)
Solve the PBE with MG.
Definition: routines.c:1341
VPUBLIC void Vclist_dtor(Vclist **thee)
Destroy object.
Definition: vclist.c:397
VPUBLIC void Vpmg_setPart(Vpmg *thee, double lowerCorner[3], double upperCorner[3], int bflags[6])
Set partition information which restricts the calculation of observables to a (rectangular) subset of...
Definition: vpmg.c:625
VPUBLIC Vpbe * Vpbe_ctor(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)
Construct Vpbe object.
Definition: vpbe.c:246
int method
Definition: mgparm.h:192
Vpmgp * pmgp
Definition: vpmg.h:119
VPUBLIC int writematMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out operator matrix from MG calculation to file.
Definition: routines.c:1653
int nion
Definition: pbeparm.h:138
VPUBLIC int setPartMG(NOsh *nosh, MGparm *mgparm, Vpmg *pmg)
Set MG partitions for calculating observables and performing I/O.
Definition: routines.c:1377
Header file for front end auxiliary routines.
VPUBLIC int Vpmg_dbForce(Vpmg *thee, double *dbForce, int atomID, Vsurf_Meth srfm)
Calculate the dielectric boundary forces on the specified atom in units of k_B T/AA.
Definition: vpmg.c:6000
VPUBLIC int Vacc_wcaEnergy(Vacc *acc, APOLparm *apolparm, Valist *alist, Vclist *clist)
Return the WCA integral energy.
Definition: vacc.c:1733
VPUBLIC void killChargeMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded charge maps.
Definition: routines.c:849
#define APBS_TIMER_SETUP
APBS setup timer ID.
Definition: vhal.h:336
Definition: vhal.h:211
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
VPUBLIC void killForce(Vmem *mem, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Free memory from MG force calculation.
Definition: routines.c:1636
Definition: pbeparm.h:82
double savForce[3]
Definition: routines.h:88
APOLparm_calcForce calcforce
Definition: apolparm.h:170
double sdie
Definition: pbeparm.h:148
double temp
Definition: pbeparm.h:156
char dielXpath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:228
VPUBLIC int loadKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the kappa maps given in NOsh into grid objects.
Definition: routines.c:575
VPUBLIC double Vacc_totalSASA(Vacc *thee, double radius)
Return the total solvent accessible surface area (SASA)
Definition: vacc.c:774
VPUBLIC int loadDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Load the dielectric maps given in NOsh into grid objects.
Definition: routines.c:250
double center[3]
Definition: mgparm.h:138
int writematflag
Definition: pbeparm.h:196
VPUBLIC int Vacc_wcaForceAtom(Vacc *thee, APOLparm *apolparm, Vclist *clist, Vatom *atom, double *force)
Return the WCA integral force.
Definition: vacc.c:1768
int ncalc
Definition: nosh.h:188
VPUBLIC Vrc_Codes Vfetk_loadMesh(Vfetk *thee, double center[3], double length[3], Vfetk_MeshLoad meshType, Vio *sock)
Loads a mesh into the Vfetk (and associated) object(s).
Definition: vfetk.c:980
VPUBLIC int initMG(int icalc, NOsh *nosh, MGparm *mgparm, PBEparm *pbeparm, double realCenter[3], Vpbe *pbe[NOSH_MAXCALC], Valist *alist[NOSH_MAXMOL], Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL], Vgrid *kappaMap[NOSH_MAXMOL], Vgrid *chargeMap[NOSH_MAXMOL], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC], Vgrid *potMap[NOSH_MAXMOL])
Initialize an MG calculation.
Definition: routines.c:1074
double press
Definition: apolparm.h:148
#define APBS_TIMER_SOLVER
APBS solver timer ID.
Definition: vhal.h:342
double sdens
Definition: pbeparm.h:146
Vsurf_Meth srfm
Definition: pbeparm.h:150
int useChargeMap
Definition: pbeparm.h:131
Definition: vfetk.h:123
double etol
Definition: femparm.h:142
VPUBLIC void storeAtomEnergy(Vpmg *pmg, int icalc, double **atomEnergy, int *nenergy)
Store energy in arrays for future use.
Definition: routines.c:1724
double glen[3]
Definition: mgparm.h:135
VPUBLIC void killMG(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vpmgp *pmgp[NOSH_MAXCALC], Vpmg *pmg[NOSH_MAXCALC])
Kill structures initialized during an MG calculation.
Definition: routines.c:1315
AtomData sub-class; stores atom data.
Definition: vparam.h:92
VPUBLIC int printApolForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
Definition: routines.c:3314
Definition: vhal.h:315
Definition: pbeparm.h:98
APOLparm * apolparm
Definition: nosh.h:168
VPUBLIC void printFEPARM(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Print out FE-specific params loaded from input.
Definition: routines.c:3743
VPUBLIC void Vfetk_dtor(Vfetk **thee)
Object destructor.
Definition: vfetk.c:625
VPUBLIC double Vpmg_qfAtomEnergy(Vpmg *thee, Vatom *atom)
Get the per-atom "fixed charge" contribution to the electrostatic energy.
Definition: vpmg.c:1796
VPUBLIC void killPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded potential maps.
Definition: routines.c:751
VPUBLIC void Vfetk_setParameters(Vfetk *thee, PBEparm *pbeparm, FEMparm *feparm)
Set the parameter objects.
Definition: vfetk.c:615
int nkappa
Definition: nosh.h:235
VPUBLIC void killFE(NOsh *nosh, Vpbe *pbe[NOSH_MAXCALC], Vfetk *fetk[NOSH_MAXCALC], Gem *gm[NOSH_MAXMOL])
Kill structures initialized during an FE calculation.
Definition: routines.c:3526
Vdata_Format potfmt[NOSH_MAXMOL]
Definition: nosh.h:240
double ionr[MAXION]
Definition: pbeparm.h:142
double zmem
Definition: pbeparm.h:174
PBEparm * pbeparm
Definition: nosh.h:167
VPUBLIC Vrc_Codes Valist_readXML(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from an XML file.
Definition: valist.c:714
VPUBLIC void Vpmg_dtor(Vpmg **thee)
Object destructor.
Definition: vpmg.c:559
Definition: vhal.h:317
MGparm_CalcType type
Definition: mgparm.h:116
int npot
Definition: nosh.h:238
VPUBLIC int printElecEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data.
Definition: routines.c:2696
double dpos
Definition: apolparm.h:145
int proc_rank
Definition: nosh.h:203
Vdata_Format kappafmt[NOSH_MAXMOL]
Definition: nosh.h:237
Definition: vhal.h:143
Definition: vhal.h:279
double dbForce[3]
Definition: routines.h:86
VPUBLIC int partFE(int icalc, NOsh *nosh, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Partition mesh (if applicable)
Definition: routines.c:3879
VPUBLIC void killKappaMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Destroy the loaded kappa maps.
Definition: routines.c:662
Definition: vhal.h:281
APOLparm_calcEnergy calcenergy
Definition: apolparm.h:167
#define VEMBED(rctag)
Allows embedding of RCS ID tags in object files.
Definition: vhal.h:563
int nelec
Definition: nosh.h:193
VPUBLIC double Vacc_totalSAV(Vacc *thee, Vclist *clist, APOLparm *apolparm, double radius)
Return the total solvent accessible volume (SAV)
Definition: vacc.c:1515
enum eMGparm_CalcType MGparm_CalcType
Declare MGparm_CalcType type.
Definition: mgparm.h:89
int ncharge
Definition: nosh.h:241
double mdie
Definition: pbeparm.h:178
VPUBLIC void Vacc_atomdSASA(Vacc *thee, double dpos, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible area.
Definition: vacc.c:1332
Definition: vhal.h:314
Definition: vhal.h:312
int dielMapID
Definition: pbeparm.h:123
VPUBLIC int Vparam_readFlatFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read a flat-file format parameter database.
Definition: vparam.c:445
VPUBLIC int Vpmg_fillArray(Vpmg *thee, double *vec, Vdata_Type type, double parm, Vhal_PBEType pbetype, PBEparm *pbeparm)
Fill the specified array with accessibility values.
Definition: vpmg.c:890
VPUBLIC int Vstring_strcasecmp(const char *s1, const char *s2)
Case-insensitive string comparison (BSD standard)
Definition: vstring.c:66
double * rwork
Definition: vpmg.h:137
VPUBLIC int energyFE(NOsh *nosh, int icalc, Vfetk *fetk[NOSH_MAXCALC], int *nenergy, double *totEnergy, double *qfEnergy, double *qmEnergy, double *dielEnergy)
Calculate electrostatic energies from FE solution.
Definition: routines.c:4019
VPUBLIC Vatom * Valist_getAtom(Valist *thee, int i)
Get pointer to particular atom in list.
Definition: valist.c:115
VPUBLIC int printForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data (deprecated...see printElecForce)
Definition: routines.c:2823
int ndiel
Definition: nosh.h:227
VPUBLIC void printPBEPARM(PBEparm *pbeparm)
Print out generic PBE params loaded from input.
Definition: routines.c:867
#define Vunit_Na
Avogadro&#39;s number.
Definition: vunit.h:100
VPUBLIC int writedataMG(int rank, NOsh *nosh, PBEparm *pbeparm, Vpmg *pmg)
Write out observables from MG calculation to file.
Definition: routines.c:2237
VPUBLIC Vparam_AtomData * Vparam_getAtomData(Vparam *thee, char resName[VMAX_ARGLEN], char atomName[VMAX_ARGLEN])
Get atom data.
Definition: vparam.c:267
Parameter structure for FEM-specific variables from input files.
Definition: femparm.h:133
double grid[3]
Definition: mgparm.h:133
Definition: nosh.h:130
double smsize
Definition: pbeparm.h:159
int nonlintype
Definition: mgparm.h:189
Structure to hold atomic forces.
Definition: routines.h:83
Reads and assigns charge/radii parameters.
Definition: vparam.h:135
double gamma
Definition: apolparm.h:163
int maxsolve
Definition: femparm.h:165
Definition: vfetk.h:124
Definition: vhal.h:277
Definition: vhal.h:218
int printcalc[NOSH_MAXPRINT][NOSH_MAXPOP]
Definition: nosh.h:251
VPUBLIC double Vpmg_dielEnergy(Vpmg *thee, int extFlag)
Get the "polarization" contribution to the electrostatic energy.
Definition: vpmg.c:1284
VPUBLIC Vpmg * Vpmg_ctor(Vpmgp *pmgp, Vpbe *pbe, int focusFlag, Vpmg *pmgOLD, MGparm *mgparm, PBEparm_calcEnergy energyFlag)
Constructor for the Vpmg class (allocates new memory)
Definition: vpmg.c:141
Contains public data members for Vpmgp class/module.
Definition: vpmgp.h:80
#define Vunit_kb
Boltzmann constant.
Definition: vunit.h:96
PBEparm_calcEnergy calcenergy
Definition: pbeparm.h:165
Definition: vhal.h:275
VPUBLIC int printEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Combine and pretty-print energy data (deprecated...see printElecEnergy)
Definition: routines.c:2628
int setwat
Definition: apolparm.h:180
VPUBLIC double returnEnergy(Vcom *com, NOsh *nosh, double totEnergy[NOSH_MAXCALC], int iprint)
Access net local energy.
Definition: routines.c:2596
#define APBS_TIMER_ENERGY
APBS energy timer ID.
Definition: vhal.h:348
VPUBLIC int loadPotMaps(NOsh *nosh, Vgrid *map[NOSH_MAXMOL])
Load the potential maps given in NOsh into grid objects.
Definition: routines.c:679
Definition: vhal.h:313
NOsh_ParmFormat parmfmt
Definition: nosh.h:226
Definition: vfetk.h:90
double ionc[MAXION]
Definition: pbeparm.h:141
int pdime[3]
Definition: mgparm.h:178
NOsh_MolFormat molfmt[NOSH_MAXMOL]
Definition: nosh.h:221
char chargepath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:242
Parameter structure for PBE variables from input files.
Definition: pbeparm.h:117
int ny
Definition: vpmgp.h:84
Definition: nosh.h:131
Definition: vhal.h:213
VPUBLIC int Vpmg_fillco(Vpmg *thee, Vsurf_Meth surfMeth, double splineWin, Vchrg_Meth chargeMeth, int useDielXMap, Vgrid *dielXMap, int useDielYMap, Vgrid *dielYMap, int useDielZMap, Vgrid *dielZMap, int useKappaMap, Vgrid *kappaMap, int usePotMap, Vgrid *potMap, int useChargeMap, Vgrid *chargeMap)
Fill the coefficient arrays prior to solving the equation.
Definition: vpmg.c:5645
Definition: vfetk.h:158
Definition: vhal.h:142
double sasa
Definition: apolparm.h:175
int nx
Definition: vpmgp.h:83
VPUBLIC void killMolecules(NOsh *nosh, Valist *alist[NOSH_MAXMOL])
Destroy the loaded molecules.
Definition: routines.c:233
double watepsilon
Definition: apolparm.h:174
PBEparm_calcForce calcforce
Definition: pbeparm.h:167
Vbcfl bcfl
Definition: pbeparm.h:136
int printnarg[NOSH_MAXPRINT]
Definition: nosh.h:250
VPUBLIC int forceMG(Vmem *mem, NOsh *nosh, PBEparm *pbeparm, MGparm *mgparm, Vpmg *pmg, int *nforce, AtomForce **atomForce, Valist *alist[NOSH_MAXMOL])
Calculate forces from MG solution.
Definition: routines.c:1490
double pdie
Definition: pbeparm.h:144
double partDisjLength[3]
Definition: mgparm.h:173
Surface object list of per-atom surface points.
Definition: vacc.h:84
VPUBLIC Vrc_Codes Valist_readPDB(Valist *thee, Vparam *param, Vio *sock)
Fill atom list with information from a PDB file.
Definition: valist.c:515
char dielZpath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:232
FEMparm_EstType akeySOLVE
Definition: femparm.h:152
Parameter structure for MG-specific variables from input files.
Definition: mgparm.h:114
Vdata_Format chargefmt[NOSH_MAXMOL]
Definition: nosh.h:243
int nmol
Definition: nosh.h:219
VPUBLIC void Vpbe_dtor(Vpbe **thee)
Object destructor.
Definition: vpbe.c:467
FEMparm_EtolType ekey
Definition: femparm.h:145
int bogus
Definition: nosh.h:205
VPUBLIC void killDielMaps(NOsh *nosh, Vgrid *dielXMap[NOSH_MAXMOL], Vgrid *dielYMap[NOSH_MAXMOL], Vgrid *dielZMap[NOSH_MAXMOL])
Destroy the loaded dielectric.
Definition: routines.c:550
VPUBLIC Vrc_Codes initFE(int icalc, NOsh *nosh, FEMparm *feparm, PBEparm *pbeparm, Vpbe *pbe[NOSH_MAXCALC], Valist *alist[NOSH_MAXMOL], Vfetk *fetk[NOSH_MAXCALC])
Initialize FE solver objects.
Definition: routines.c:3554
double ionq[MAXION]
Definition: pbeparm.h:140
Definition: vfetk.h:87
AM * am
Definition: vfetk.h:182
VPUBLIC int printElecForce(Vcom *com, NOsh *nosh, int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC], int iprint)
Combine and pretty-print force data.
Definition: routines.c:3071
Definition: vhal.h:283
int useDielMap
Definition: pbeparm.h:121
double memv
Definition: pbeparm.h:180
int useMesh
Definition: femparm.h:173
double ofrac
Definition: mgparm.h:184
VPUBLIC double * Vatom_getPosition(Vatom *thee)
Get atomic position.
Definition: vatom.c:63
VPUBLIC void printMGPARM(MGparm *mgparm, double realCenter[3])
Print out MG-specific params loaded from input.
Definition: routines.c:1041
int dime[3]
Definition: mgparm.h:120
VPUBLIC int printApolEnergy(NOsh *nosh, int iprint)
Combine and pretty-print energy data.
Definition: routines.c:2761
int kappaMapID
Definition: pbeparm.h:126
Contains public data members for Vatom class/module.
Definition: vatom.h:84
char writematstem[VMAX_ARGLEN]
Definition: pbeparm.h:195
VPUBLIC int Vgrid_readDX(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read in data in OpenDX grid format.
Definition: vgrid.c:575
Container class for list of atom objects.
Definition: valist.h:78
Valist * alist
Definition: vpbe.h:88
VPUBLIC Vpmgp * Vpmgp_ctor(MGparm *mgparm)
Construct PMG parameter object and initialize to default values.
Definition: vpmgp.c:76
Vchrg_Meth chgm
Definition: mgparm.h:122
double ycent
Definition: vpmgp.h:119
Class for parsing fixed format input files.
Definition: nosh.h:183
char meshpath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:245
int meth
Definition: vpmgp.h:144
double srad
Definition: pbeparm.h:152
double sasaForce[3]
Definition: routines.h:87
VPUBLIC int Vparam_readXMLFile(Vparam *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname)
Read an XML format parameter database.
Definition: vparam.c:306
int elec2calc[NOSH_MAXCALC]
Definition: nosh.h:209
Contains public data members for Vfetk class/module.
Definition: vfetk.h:176
VPUBLIC int Vpmg_ibForce(Vpmg *thee, double *force, int atomID, Vsurf_Meth srfm)
Calculate the osmotic pressure on the specified atom in units of k_B T/AA.
Definition: vpmg.c:5835
VPUBLIC int loadMolecules(NOsh *nosh, Vparam *param, Valist *alist[NOSH_MAXMOL])
Load the molecules given in NOsh into atom lists.
Definition: routines.c:95
VPUBLIC void Vgrid_writeGZ(Vgrid *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, char *title, double *pvec)
Write out OpenDX data in GZIP format.
Definition: vgrid.c:807
int number
Definition: valist.h:80
Definition: vfetk.h:88
double watsigma
Definition: apolparm.h:173
int writemat
Definition: pbeparm.h:191
#define NOSH_MAXCALC
Maximum number of calculations in a run.
Definition: nosh.h:83
int partDisjOwnSide[6]
Definition: mgparm.h:175
Vpbe * pbe
Definition: vpmg.h:120
VPUBLIC int forceAPOL(Vacc *acc, Vmem *mem, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist, Vclist *clist)
Calculate non-polar forces.
Definition: routines.c:4562
Vdata_Type writetype[PBEPARM_MAXWRITE]
Definition: pbeparm.h:188
VPUBLIC int Vpmg_qfForce(Vpmg *thee, double *force, int atomID, Vchrg_Meth chgm)
Calculate the "charge-field" force on the specified atom in units of k_B T/AA.
Definition: vpmg.c:6257
VPUBLIC int solveFE(int icalc, PBEparm *pbeparm, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Solve-estimate-refine.
Definition: routines.c:3956
VPUBLIC Vrc_Codes Valist_readPQR(Valist *thee, Vparam *params, Vio *sock)
Fill atom list with information from a PQR file.
Definition: valist.c:606
VPUBLIC double Vfetk_energy(Vfetk *thee, int color, int nonlin)
Return the total electrostatic energy.
Definition: vfetk.c:693
int numwrite
Definition: pbeparm.h:185
Atom cell list.
Definition: vclist.h:117
enum eVfetk_MeshLoad Vfetk_MeshLoad
Declare FEMparm_GuessType type.
Definition: vfetk.h:114
VPUBLIC void startVio()
Wrapper to start MALOC Vio layer.
Definition: routines.c:58
int maxvert
Definition: femparm.h:167
double wcaForce[3]
Definition: routines.h:89
int useAqua
Definition: mgparm.h:195
double partDisjCenter[3]
Definition: mgparm.h:171
char parmpath[VMAX_ARGLEN]
Definition: nosh.h:225
double qfForce[3]
Definition: routines.h:85
VPUBLIC int Vpmg_solve(Vpmg *thee)
Solve the PBE using PMG.
Definition: vpmg.c:399
int ispara
Definition: nosh.h:202
int nlev
Definition: mgparm.h:128
VPUBLIC double Vpmg_energy(Vpmg *thee, int extFlag)
Get the total electrostatic energy.
Definition: vpmg.c:1253
VPUBLIC double Vpbe_getDeblen(Vpbe *thee)
Get Debye-Huckel screening length.
Definition: vpbe.c:141
int printop[NOSH_MAXPRINT][NOSH_MAXPOP]
Definition: nosh.h:252
char molpath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:220
VPUBLIC int Vfetk_fillArray(Vfetk *thee, Bvec *vec, Vdata_Type type)
Fill an array with the specified data.
Definition: vfetk.c:2299
int usePotMap
Definition: pbeparm.h:127
VPUBLIC double Vatom_getRadius(Vatom *thee)
Get atomic position.
Definition: vatom.c:105
Definition: nosh.h:100
Vmem * mem
Definition: vacc.h:110
Valist * alist
Definition: vacc.h:111
VPUBLIC double Vatom_getCharge(Vatom *thee)
Get atomic charge.
Definition: vatom.c:119
char writestem[PBEPARM_MAXWRITE][VMAX_ARGLEN]
Definition: pbeparm.h:186
double wcaEnergy
Definition: apolparm.h:177
VPUBLIC void Valist_dtor(Valist **thee)
Destroys atom list object.
Definition: valist.c:167
FEMparm * femparm
Definition: nosh.h:166
char apolname[NOSH_MAXCALC][VMAX_ARGLEN]
Definition: nosh.h:257
VPUBLIC void Vacc_atomdSAV(Vacc *thee, double srad, Vatom *atom, double *dSA)
Get the derivatve of solvent accessible volume.
Definition: vacc.c:1212
VPUBLIC int initAPOL(NOsh *nosh, Vmem *mem, Vparam *param, APOLparm *apolparm, int *nforce, AtomForce **atomForce, Valist *alist)
Upperlevel routine to the non-polar energy and force routines.
Definition: routines.c:4308
double zmin
Definition: vpmgp.h:194
VPUBLIC int Vfetk_write(Vfetk *thee, const char *iodev, const char *iofmt, const char *thost, const char *fname, Bvec *vec, Vdata_Format format)
Write out data.
Definition: vfetk.c:2464
VPUBLIC Valist * Valist_ctor()
Construct the atom list object.
Definition: valist.c:138
#define APBS_TIMER_FORCE
APBS force timer ID.
Definition: vhal.h:354
double maxIonRadius
Definition: vpbe.h:100
Parameter structure for APOL-specific variables from input files.
Definition: apolparm.h:129
VPUBLIC int writedataXML(NOsh *nosh, Vcom *com, const char *fname, double totEnergy[NOSH_MAXCALC], double qfEnergy[NOSH_MAXCALC], double qmEnergy[NOSH_MAXCALC], double dielEnergy[NOSH_MAXCALC], int nenergy[NOSH_MAXCALC], double *atomEnergy[NOSH_MAXCALC], int nforce[NOSH_MAXCALC], AtomForce *atomForce[NOSH_MAXCALC])
Write out information to an XML file.
Definition: routines.c:1977
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:873
int nmesh
Definition: nosh.h:244
VPUBLIC int energyMG(NOsh *nosh, int icalc, Vpmg *pmg, int *nenergy, double *totEnergy, double *qfEnergy, double *qmEnergy, double *dielEnergy)
Calculate electrostatic energies from MG solution.
Definition: routines.c:1423
FEMparm_EstType akeyPRE
Definition: femparm.h:148
double * pvec
Definition: vpmg.h:154
VPUBLIC int preRefineFE(int icalc, FEMparm *feparm, Vfetk *fetk[NOSH_MAXCALC])
Pre-refine mesh before solve.
Definition: routines.c:3886
double xmax
Definition: vpmgp.h:195
VPUBLIC int writedataFE(int rank, NOsh *nosh, PBEparm *pbeparm, Vfetk *fetk)
Write FEM data to files.
Definition: routines.c:4140
VPUBLIC Vgrid * Vgrid_ctor(int nx, int ny, int nz, double hx, double hy, double hzed, double xmin, double ymin, double zmin, double *data)
Construct Vgrid object with values obtained from Vpmg_readDX (for example)
Definition: vgrid.c:77
VPUBLIC void killEnergy()
Kill arrays allocated for energies.
Definition: routines.c:1628
double radius
Definition: vparam.h:96
int nz
Definition: vpmgp.h:85
char dielYpath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:230
char kappapath[NOSH_MAXMOL][VMAX_ARGLEN]
Definition: nosh.h:236