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