Bug Summary

File:gromacs/gmxana/gmx_dos.c
Location:line 110, column 9
Description:Value stored to 'fd0' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
8 *
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
13 *
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
23 *
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
31 *
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
34 */
35#ifdef HAVE_CONFIG_H1
36#include <config.h>
37#endif
38
39#include <math.h>
40#include <stdio.h>
41#include <stdlib.h>
42#include <string.h>
43
44#include "gromacs/fileio/confio.h"
45#include "copyrite.h"
46#include "gromacs/utility/fatalerror.h"
47#include "gromacs/utility/futil.h"
48#include "gstat.h"
49#include "macros.h"
50#include "gromacs/math/utilities.h"
51#include "physics.h"
52#include "index.h"
53#include "gromacs/utility/smalloc.h"
54#include "gromacs/commandline/pargs.h"
55#include "txtdump.h"
56#include "typedefs.h"
57#include "gromacs/math/vec.h"
58#include "gromacs/fileio/xvgr.h"
59#include "viewit.h"
60#include "correl.h"
61#include "gmx_ana.h"
62#include "gromacs/fft/fft.h"
63#include "gromacs/fileio/trxio.h"
64
65enum {
66 VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
67};
68
69static double FD(double Delta, double f)
70{
71 return (2*pow(Delta, -4.5)*pow(f, 7.5) -
72 6*pow(Delta, -3)*pow(f, 5) -
73 pow(Delta, -1.5)*pow(f, 3.5) +
74 6*pow(Delta, -1.5)*pow(f, 2.5) +
75 2*f - 2);
76}
77
78static double YYY(double f, double y)
79{
80 return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
81 (2+6*y)*f - 2);
82}
83
84static double calc_compress(double y)
85{
86 if (y == 1)
87 {
88 return 0;
89 }
90 return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
91}
92
93static double bisector(double Delta, double tol,
94 double ff0, double ff1,
95 double ff(double, double))
96{
97 double fd0, fd, fd1, f, f0, f1;
98 double tolmin = 1e-8;
99
100 f0 = ff0;
101 f1 = ff1;
102 if (tol < tolmin)
103 {
104 fprintf(stderrstderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
105 tol = tolmin;
106 }
107
108 do
109 {
110 fd0 = ff(Delta, f0);
Value stored to 'fd0' is never read
111 fd1 = ff(Delta, f1);
112 f = (f0+f1)*0.5;
113 fd = ff(Delta, f);
114 if (fd < 0)
115 {
116 f0 = f;
117 }
118 else if (fd > 0)
119 {
120 f1 = f;
121 }
122 else
123 {
124 return f;
125 }
126 }
127 while ((f1-f0) > tol);
128
129 return f;
130}
131
132static double calc_fluidicity(double Delta, double tol)
133{
134 return bisector(Delta, tol, 0, 1, FD);
135}
136
137static double calc_y(double f, double Delta, double toler)
138{
139 double y1, y2;
140
141 y1 = pow(f/Delta, 1.5);
142 y2 = bisector(f, toler, 0, 10000, YYY);
143 if (fabs((y1-y2)/(y1+y2)) > 100*toler)
144 {
145 fprintf(stderrstderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
146 y1, y2);
147 }
148
149 return y1;
150}
151
152static double calc_Shs(double f, double y)
153{
154 double fy = f*y;
155
156 return BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
157}
158
159static real wCsolid(real nu, real beta)
160{
161 real bhn = beta*PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu;
162 real ebn, koko;
163
164 if (bhn == 0)
165 {
166 return 1.0;
167 }
168 else
169 {
170 ebn = exp(bhn);
171 koko = sqr(1-ebn);
172 return sqr(bhn)*ebn/koko;
173 }
174}
175
176static real wSsolid(real nu, real beta)
177{
178 real bhn = beta*PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu;
179
180 if (bhn == 0)
181 {
182 return 1;
183 }
184 else
185 {
186 return bhn/(exp(bhn)-1) - log(1-exp(-bhn));
187 }
188}
189
190static real wAsolid(real nu, real beta)
191{
192 real bhn = beta*PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu;
193
194 if (bhn == 0)
195 {
196 return 0;
197 }
198 else
199 {
200 return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
201 }
202}
203
204static real wEsolid(real nu, real beta)
205{
206 real bhn = beta*PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu;
207
208 if (bhn == 0)
209 {
210 return 1;
211 }
212 else
213 {
214 return bhn/2 + bhn/(exp(bhn)-1)-1;
215 }
216}
217
218static void dump_fy(output_env_t oenv, real toler)
219{
220 FILE *fp;
221 double Delta, f, y, DD;
222 const char *leg[] = { "f", "fy", "y" };
223
224 DD = pow(10.0, 0.125);
225 fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);
226 xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv);
227 fprintf(fp, "@ world 1e-05, 0, 1000, 1\n");
228 fprintf(fp, "@ xaxes scale Logarithmic\n");
229 for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
230 {
231 f = calc_fluidicity(Delta, toler);
232 y = calc_y(f, Delta, toler);
233 fprintf(fp, "%10g %10g %10g %10g\n", Delta, f, f*y, y);
234 }
235 xvgrclose(fp);
236}
237
238static void dump_w(output_env_t oenv, real beta)
239{
240 FILE *fp;
241 double nu;
242 const char *leg[] = { "wCv", "wS", "wA", "wE" };
243
244 fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
245 "w", oenv);
246 xvgr_legend(fp, asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv);
247 for (nu = 1; (nu < 100); nu += 0.05)
248 {
249 fprintf(fp, "%10g %10g %10g %10g %10g\n", beta*PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu,
250 wCsolid(nu, beta), wSsolid(nu, beta),
251 wAsolid(nu, beta), wEsolid(nu, beta));
252 }
253 xvgrclose(fp);
254}
255
256int gmx_dos(int argc, char *argv[])
257{
258 const char *desc[] = {
259 "[THISMODULE] computes the Density of States from a simulations.",
260 "In order for this to be meaningful the velocities must be saved",
261 "in the trajecotry with sufficiently high frequency such as to cover",
262 "all vibrations. For flexible systems that would be around a few fs",
263 "between saving. Properties based on the DoS are printed on the",
264 "standard output."
265 };
266 const char *bugs[] = {
267 "This program needs a lot of memory: total usage equals the number of atoms times 3 times number of frames times 4 (or 8 when run in double precision)."
268 };
269 FILE *fp, *fplog;
270 t_topology top;
271 int ePBC = -1;
272 t_trxframe fr;
273 matrix box;
274 int gnx;
275 char title[256];
276 real t0, t1, m;
277 t_trxstatus *status;
278 int nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
279 double rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
280 real **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
281 output_env_t oenv;
282 gmx_fft_t fft;
283 double cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
284 double wCdiff, wSdiff, wAdiff, wEdiff;
285
286 static gmx_bool bVerbose = TRUE1, bAbsolute = FALSE0, bNormalize = FALSE0;
287 static gmx_bool bRecip = FALSE0, bDump = FALSE0;
288 static real Temp = 298.15, toler = 1e-6;
289 t_pargs pa[] = {
290 { "-v", FALSE0, etBOOL, {&bVerbose},
291 "Be loud and noisy." },
292 { "-recip", FALSE0, etBOOL, {&bRecip},
293 "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
294 { "-abs", FALSE0, etBOOL, {&bAbsolute},
295 "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
296 { "-normdos", FALSE0, etBOOL, {&bNormalize},
297 "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
298 { "-T", FALSE0, etREAL, {&Temp},
299 "Temperature in the simulation" },
300 { "-toler", FALSE0, etREAL, {&toler},
301 "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
302 { "-dump", FALSE0, etBOOL, {&bDump},
303 "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
304 };
305
306 t_filenm fnm[] = {
307 { efTRN, "-f", NULL((void*)0), ffREAD1<<1 },
308 { efTPX, "-s", NULL((void*)0), ffREAD1<<1 },
309 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
310 { efXVG, "-vacf", "vacf", ffWRITE1<<2 },
311 { efXVG, "-mvacf", "mvacf", ffWRITE1<<2 },
312 { efXVG, "-dos", "dos", ffWRITE1<<2 },
313 { efLOG, "-g", "dos", ffWRITE1<<2 },
314 };
315#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
316 int npargs;
317 t_pargs *ppa;
318 const char *DoSlegend[] = {
319 "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
320 };
321
322 npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])));
323 ppa = add_acf_pargs(&npargs, pa);
324 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_BE_NICE(1<<13),
325 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, ppa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc,
326 asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs, &oenv))
327 {
328 return 0;
329 }
330
331 beta = 1/(Temp*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3)));
332 if (bDump)
333 {
334 printf("Dumping reference figures. Thanks for your patience.\n");
335 dump_fy(oenv, toler);
336 dump_w(oenv, beta);
337 exit(0);
338 }
339
340 fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w");
341 fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
342 please_cite(fplog, "Pascal2011a");
343 please_cite(fplog, "Caleman2011b");
344
345 read_tps_conf(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, NULL((void*)0), NULL((void*)0), box,
346 TRUE1);
347 V = det(box);
348 tmass = 0;
349 for (i = 0; (i < top.atoms.nr); i++)
350 {
351 tmass += top.atoms.atom[i].m;
352 }
353
354 Natom = top.atoms.nr;
355 Nmol = top.mols.nr;
356 gnx = Natom*DIM3;
357
358 /* Correlation stuff */
359 snew(c1, gnx)(c1) = save_calloc("c1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 359, (gnx), sizeof(*(c1)));
360 for (i = 0; (i < gnx); i++)
361 {
362 c1[i] = NULL((void*)0);
363 }
364
365 read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &fr, TRX_NEED_V(1<<3));
366 t0 = fr.time;
367
368 n_alloc = 0;
369 nframes = 0;
370 Vsum = V2sum = 0;
371 nV = 0;
372 do
373 {
374 if (fr.bBox)
375 {
376 V = det(fr.box);
377 V2sum += V*V;
378 Vsum += V;
379 nV++;
380 }
381 if (nframes >= n_alloc)
382 {
383 n_alloc += 100;
384 for (i = 0; i < gnx; i++)
385 {
386 srenew(c1[i], n_alloc)(c1[i]) = save_realloc("c1[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 386, (c1[i]), (n_alloc), sizeof(*(c1[i])));
387 }
388 }
389 for (i = 0; i < gnx; i += DIM3)
390 {
391 c1[i+XX0][nframes] = fr.v[i/DIM3][XX0];
392 c1[i+YY1][nframes] = fr.v[i/DIM3][YY1];
393 c1[i+ZZ2][nframes] = fr.v[i/DIM3][ZZ2];
394 }
395
396 t1 = fr.time;
397
398 nframes++;
399 }
400 while (read_next_frame(oenv, status, &fr));
401
402 close_trj(status);
403
404 dt = (t1-t0)/(nframes-1);
405 if (nV > 0)
406 {
407 V = Vsum/nV;
408 }
409 if (bVerbose)
410 {
411 printf("Going to do %d fourier transforms of length %d. Hang on.\n",
412 gnx, nframes);
413 }
414 low_do_autocorr(NULL((void*)0), oenv, NULL((void*)0), nframes, gnx, nframes, c1, dt, eacNormal(1<<0), 0, FALSE0,
415 FALSE0, FALSE0, -1, -1, 0);
416 snew(dos, DOS_NR)(dos) = save_calloc("dos", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 416, (DOS_NR), sizeof(*(dos)));
417 for (j = 0; (j < DOS_NR); j++)
418 {
419 snew(dos[j], nframes+4)(dos[j]) = save_calloc("dos[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 419, (nframes+4), sizeof(*(dos[j])));
420 }
421
422 if (bVerbose)
423 {
424 printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
425 }
426 for (i = 0; (i < gnx); i += DIM3)
427 {
428 mi = top.atoms.atom[i/DIM3].m;
429 for (j = 0; (j < nframes/2); j++)
430 {
431 c1j = (c1[i+XX0][j] + c1[i+YY1][j] + c1[i+ZZ2][j]);
432 dos[VACF][j] += c1j/Natom;
433 dos[MVACF][j] += mi*c1j;
434 }
435 }
436 fp = xvgropen(opt2fn("-vacf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Velocity ACF",
437 "Time (ps)", "C(t)", oenv);
438 snew(tt, nframes/2)(tt) = save_calloc("tt", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 438, (nframes/2), sizeof(*(tt)));
439 for (j = 0; (j < nframes/2); j++)
440 {
441 tt[j] = j*dt;
442 fprintf(fp, "%10g %10g\n", tt[j], dos[VACF][j]);
443 }
444 xvgrclose(fp);
445 fp = xvgropen(opt2fn("-mvacf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Mass-weighted velocity ACF",
446 "Time (ps)", "C(t)", oenv);
447 for (j = 0; (j < nframes/2); j++)
448 {
449 fprintf(fp, "%10g %10g\n", tt[j], dos[MVACF][j]);
450 }
451 xvgrclose(fp);
452
453 if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
454 GMX_FFT_FLAG_NONE)) != 0)
455 {
456 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 456, "gmx_fft_init_1d_real returned %d", fftcode);
457 }
458 if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
459 (void *)dos[MVACF], (void *)dos[DOS])) != 0)
460 {
461 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 461, "gmx_fft_1d_real returned %d", fftcode);
462 }
463
464 /* First compute the DoS */
465 /* Magic factor of 8 included now. */
466 bfac = 8*dt*beta/2;
467 dos2 = 0;
468 snew(nu, nframes/4)(nu) = save_calloc("nu", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dos.c"
, 468, (nframes/4), sizeof(*(nu)));
469 for (j = 0; (j < nframes/4); j++)
470 {
471 nu[j] = 2*j/(t1-t0);
472 dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
473 if (bAbsolute)
474 {
475 dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
476 }
477 else
478 {
479 dos[DOS][j] = bfac*dos[DOS][2*j];
480 }
481 }
482 /* Normalize it */
483 dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL((void*)0), nframes/4, &stddev);
484 if (bNormalize)
485 {
486 for (j = 0; (j < nframes/4); j++)
487 {
488 dos[DOS][j] *= 3*Natom/dostot;
489 }
490 }
491
492 /* Now analyze it */
493 DoS0 = dos[DOS][0];
494
495 /* Note this eqn. is incorrect in Pascal2011a! */
496 Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI3.14159265358979323846*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*Temp*Natom/tmass)*
497 pow((Natom/V), 1.0/3.0)*pow(6/M_PI3.14159265358979323846, 2.0/3.0));
498 f = calc_fluidicity(Delta, toler);
499 y = calc_y(f, Delta, toler);
500 z = calc_compress(y);
501 Sig = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*(5.0/2.0+log(2*M_PI3.14159265358979323846*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*Temp/(sqr(PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))))*V/(f*Natom)));
502 Shs = Sig+calc_Shs(f, y);
503 rho = (tmass*AMU(1.6605402e-27))/(V*NANO(1e-9)*NANO(1e-9)*NANO(1e-9));
504 sigHS = pow(6*y*V/(M_PI3.14159265358979323846*Natom), 1.0/3.0);
505
506 fprintf(fplog, "System = \"%s\"\n", title);
507 fprintf(fplog, "Nmol = %d\n", Nmol);
508 fprintf(fplog, "Natom = %d\n", Natom);
509 fprintf(fplog, "dt = %g ps\n", dt);
510 fprintf(fplog, "tmass = %g amu\n", tmass);
511 fprintf(fplog, "V = %g nm^3\n", V);
512 fprintf(fplog, "rho = %g g/l\n", rho);
513 fprintf(fplog, "T = %g K\n", Temp);
514 fprintf(fplog, "beta = %g mol/kJ\n", beta);
515
516 fprintf(fplog, "\nDoS parameters\n");
517 fprintf(fplog, "Delta = %g\n", Delta);
518 fprintf(fplog, "fluidicity = %g\n", f);
519 fprintf(fplog, "hard sphere packing fraction = %g\n", y);
520 fprintf(fplog, "hard sphere compressibility = %g\n", z);
521 fprintf(fplog, "ideal gas entropy = %g\n", Sig);
522 fprintf(fplog, "hard sphere entropy = %g\n", Shs);
523 fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
524 fprintf(fplog, "DoS0 = %g\n", DoS0);
525 fprintf(fplog, "Dos2 = %g\n", dos2);
526 fprintf(fplog, "DoSTot = %g\n", dostot);
527
528 /* Now compute solid (2) and diffusive (3) components */
529 fp = xvgropen(opt2fn("-dos", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Density of states",
530 bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
531 "\\f{4}S(\\f{12}n\\f{4})", oenv);
532 xvgr_legend(fp, asize(DoSlegend)((int)(sizeof(DoSlegend)/sizeof((DoSlegend)[0]))), DoSlegend, oenv);
533 recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT(2.9979245800E05)) : 1.0;
534 for (j = 0; (j < nframes/4); j++)
535 {
536 dos[DOS_DIFF][j] = DoS0/(1+sqr(DoS0*M_PI3.14159265358979323846*nu[j]/(6*f*Natom)));
537 dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
538 fprintf(fp, "%10g %10g %10g %10g\n",
539 recip_fac*nu[j],
540 dos[DOS][j]/recip_fac,
541 dos[DOS_SOLID][j]/recip_fac,
542 dos[DOS_DIFF][j]/recip_fac);
543 }
544 xvgrclose(fp);
545
546 /* Finally analyze the results! */
547 wCdiff = 0.5;
548 wSdiff = Shs/(3*BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))); /* Is this correct? */
549 wEdiff = 0.5;
550 wAdiff = wEdiff-wSdiff;
551 for (j = 0; (j < nframes/4); j++)
552 {
553 dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
554 dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
555 dos[DOS_S][j] = (dos[DOS_DIFF][j]*wSdiff +
556 dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
557 dos[DOS_A][j] = (dos[DOS_DIFF][j]*wAdiff +
558 dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
559 dos[DOS_E][j] = (dos[DOS_DIFF][j]*wEdiff +
560 dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
561 }
562 DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL((void*)0), nframes/2, &stddev);
563 DiffCoeff = 1000*DiffCoeff/3.0;
564 fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
565 DiffCoeff);
566 fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
567 1000*DoS0/(12*tmass*beta));
568
569 cP = BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3)) * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL((void*)0),
570 nframes/4, &stddev);
571 fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
572
573 /*
574 S = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
575 nframes/4,&stddev);
576 fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
577 A = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
578 nframes/4,&stddev);
579 fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
580 E = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
581 nframes/4,&stddev);
582 fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
583 */
584 fprintf(fplog, "\nArrivederci!\n");
585 gmx_fio_fclose(fplog);
586
587 do_view(oenv, ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
588
589 return 0;
590}