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