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