File: | gromacs/gmxana/gmx_dos.c |
Location: | line 111, column 9 |
Description: | Value stored to 'fd1' is never read |
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 | |
65 | enum { |
66 | VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR |
67 | }; |
68 | |
69 | static 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 | |
78 | static 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 | |
84 | static 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 | |
93 | static 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); |
111 | fd1 = ff(Delta, f1); |
Value stored to 'fd1' is never read | |
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 | |
132 | static double calc_fluidicity(double Delta, double tol) |
133 | { |
134 | return bisector(Delta, tol, 0, 1, FD); |
135 | } |
136 | |
137 | static 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 | |
152 | static 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 | |
159 | static 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 | |
176 | static 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 | |
190 | static 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 | |
204 | static 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 | |
218 | static 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 | |
238 | static 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 | |
256 | int 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 | } |