Bug Summary

File:gromacs/gmxana/gmx_nmeig.c
Location:line 233, column 33
Description:Access to field 'ndata' results in a dereference of a null pointer (loaded from variable 'sparse_hessian')

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <string.h>
43
44#include "gromacs/commandline/pargs.h"
45#include "typedefs.h"
46#include "gromacs/utility/smalloc.h"
47#include "macros.h"
48#include "gromacs/math/vec.h"
49#include "pbc.h"
50#include "copyrite.h"
51#include "gromacs/utility/futil.h"
52#include "index.h"
53#include "mshift.h"
54#include "gromacs/fileio/xvgr.h"
55#include "gstat.h"
56#include "txtdump.h"
57#include "eigio.h"
58#include "mtop_util.h"
59#include "physics.h"
60#include "gmx_ana.h"
61
62#include "gromacs/linearalgebra/eigensolver.h"
63#include "gromacs/linearalgebra/mtxio.h"
64#include "gromacs/linearalgebra/sparsematrix.h"
65
66static double cv_corr(double nu, double T)
67{
68 double x = PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*T);
69 double ex = exp(x);
70
71 if (nu <= 0)
72 {
73 return BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*KILO(1e3);
74 }
75 else
76 {
77 return BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*KILO(1e3)*(ex*sqr(x)/sqr(ex-1) - 1);
78 }
79}
80
81static double u_corr(double nu, double T)
82{
83 double x = PLANCK(6.6262e-34*(6.0221367e23)/((1e-12)*(1e3)))*nu/(BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*T);
84 double ex = exp(x);
85
86 if (nu <= 0)
87 {
88 return BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*T;
89 }
90 else
91 {
92 return BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*T*(0.5*x - 1 + x/(ex-1));
93 }
94}
95
96static int get_nharm_mt(gmx_moltype_t *mt)
97{
98 static int harm_func[] = { F_BONDS };
99 int i, ft, nh;
100
101 nh = 0;
102 for (i = 0; (i < asize(harm_func)((int)(sizeof(harm_func)/sizeof((harm_func)[0])))); i++)
103 {
104 ft = harm_func[i];
105 nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
106 }
107 return nh;
108}
109
110static int get_nvsite_mt(gmx_moltype_t *mt)
111{
112 static int vs_func[] = {
113 F_VSITE2, F_VSITE3, F_VSITE3FD, F_VSITE3FAD,
114 F_VSITE3OUT, F_VSITE4FD, F_VSITE4FDN, F_VSITEN
115 };
116 int i, ft, nh;
117
118 nh = 0;
119 for (i = 0; (i < asize(vs_func)((int)(sizeof(vs_func)/sizeof((vs_func)[0])))); i++)
120 {
121 ft = vs_func[i];
122 nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
123 }
124 return nh;
125}
126
127static int get_nharm(gmx_mtop_t *mtop, int *nvsites)
128{
129 int j, mt, nh, nv;
130
131 nh = 0;
132 nv = 0;
133 for (j = 0; (j < mtop->nmolblock); j++)
134 {
135 mt = mtop->molblock[j].type;
136 nh += mtop->molblock[j].nmol * get_nharm_mt(&(mtop->moltype[mt]));
137 nv += mtop->molblock[j].nmol * get_nvsite_mt(&(mtop->moltype[mt]));
138 }
139 *nvsites = nv;
140 return nh;
141}
142
143static void
144nma_full_hessian(real * hess,
145 int ndim,
146 gmx_bool bM,
147 t_topology * top,
148 int begin,
149 int end,
150 real * eigenvalues,
151 real * eigenvectors)
152{
153 int i, j, k, l;
154 real mass_fac, rdum;
155 int natoms;
156
157 natoms = top->atoms.nr;
158
159 /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
160
161 if (bM)
162 {
163 for (i = 0; (i < natoms); i++)
164 {
165 for (j = 0; (j < DIM3); j++)
166 {
167 for (k = 0; (k < natoms); k++)
168 {
169 mass_fac = gmx_invsqrt(top->atoms.atom[i].m*top->atoms.atom[k].m)gmx_software_invsqrt(top->atoms.atom[i].m*top->atoms.atom
[k].m);
170 for (l = 0; (l < DIM3); l++)
171 {
172 hess[(i*DIM3+j)*ndim+k*DIM3+l] *= mass_fac;
173 }
174 }
175 }
176 }
177 }
178
179 /* call diagonalization routine. */
180
181 fprintf(stderrstderr, "\nDiagonalizing to find vectors %d through %d...\n", begin, end);
182 fflush(stderrstderr);
183
184 eigensolver(hess, ndim, begin-1, end-1, eigenvalues, eigenvectors);
185
186 /* And scale the output eigenvectors */
187 if (bM && eigenvectors != NULL((void*)0))
188 {
189 for (i = 0; i < (end-begin+1); i++)
190 {
191 for (j = 0; j < natoms; j++)
192 {
193 mass_fac = gmx_invsqrt(top->atoms.atom[j].m)gmx_software_invsqrt(top->atoms.atom[j].m);
194 for (k = 0; (k < DIM3); k++)
195 {
196 eigenvectors[i*ndim+j*DIM3+k] *= mass_fac;
197 }
198 }
199 }
200 }
201}
202
203
204
205static void
206nma_sparse_hessian(gmx_sparsematrix_t * sparse_hessian,
207 gmx_bool bM,
208 t_topology * top,
209 int neig,
210 real * eigenvalues,
211 real * eigenvectors)
212{
213 int i, j, k;
214 int row, col;
215 real mass_fac;
216 int iatom, katom;
217 int natoms;
218 int ndim;
219
220 natoms = top->atoms.nr;
221 ndim = DIM3*natoms;
222
223 /* Cannot check symmetry since we only store half matrix */
224 /* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
225
226 if (bM)
13
←
Taking true branch
→
227 {
228 for (iatom = 0; (iatom < natoms); iatom++)
14
←
Assuming 'iatom' is < 'natoms'
→
15
←
Loop condition is true. Entering loop body
→
229 {
230 for (j = 0; (j < DIM3); j++)
16
←
Loop condition is true. Entering loop body
→
231 {
232 row = DIM3*iatom+j;
233 for (k = 0; k < sparse_hessian->ndata[row]; k++)
17
←
Access to field 'ndata' results in a dereference of a null pointer (loaded from variable 'sparse_hessian')
234 {
235 col = sparse_hessian->data[row][k].col;
236 katom = col/3;
237 mass_fac = gmx_invsqrt(top->atoms.atom[iatom].m*top->atoms.atom[katom].m)gmx_software_invsqrt(top->atoms.atom[iatom].m*top->atoms
.atom[katom].m);
238 sparse_hessian->data[row][k].value *= mass_fac;
239 }
240 }
241 }
242 }
243 fprintf(stderrstderr, "\nDiagonalizing to find eigenvectors 1 through %d...\n", neig);
244 fflush(stderrstderr);
245
246 sparse_eigensolver(sparse_hessian, neig, eigenvalues, eigenvectors, 10000000);
247
248 /* Scale output eigenvectors */
249 if (bM && eigenvectors != NULL((void*)0))
250 {
251 for (i = 0; i < neig; i++)
252 {
253 for (j = 0; j < natoms; j++)
254 {
255 mass_fac = gmx_invsqrt(top->atoms.atom[j].m)gmx_software_invsqrt(top->atoms.atom[j].m);
256 for (k = 0; (k < DIM3); k++)
257 {
258 eigenvectors[i*ndim+j*DIM3+k] *= mass_fac;
259 }
260 }
261 }
262 }
263}
264
265
266
267int gmx_nmeig(int argc, char *argv[])
268{
269 const char *desc[] = {
270 "[THISMODULE] calculates the eigenvectors/values of a (Hessian) matrix,",
271 "which can be calculated with [gmx-mdrun].",
272 "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
273 "The structure is written first with t=0. The eigenvectors",
274 "are written as frames with the eigenvector number as timestamp.",
275 "The eigenvectors can be analyzed with [gmx-anaeig].",
276 "An ensemble of structures can be generated from the eigenvectors with",
277 "[gmx-nmens]. When mass weighting is used, the generated eigenvectors",
278 "will be scaled back to plain Cartesian coordinates before generating the",
279 "output. In this case, they will no longer be exactly orthogonal in the",
280 "standard Cartesian norm, but in the mass-weighted norm they would be.[PAR]",
281 "This program can be optionally used to compute quantum corrections to heat capacity",
282 "and enthalpy by providing an extra file argument [TT]-qcorr[tt]. See the GROMACS",
283 "manual, Chapter 1, for details. The result includes subtracting a harmonic",
284 "degree of freedom at the given temperature.",
285 "The total correction is printed on the terminal screen.",
286 "The recommended way of getting the corrections out is:[PAR]",
287 "[TT]gmx nmeig -s topol.tpr -f nm.mtx -first 7 -last 10000 -T 300 -qc [-constr][tt][PAR]",
288 "The [TT]-constr[tt] option should be used when bond constraints were used during the",
289 "simulation [BB]for all the covalent bonds[bb]. If this is not the case, ",
290 "you need to analyze the [TT]quant_corr.xvg[tt] file yourself.[PAR]",
291 "To make things more flexible, the program can also take virtual sites into account",
292 "when computing quantum corrections. When selecting [TT]-constr[tt] and",
293 "[TT]-qc[tt], the [TT]-begin[tt] and [TT]-end[tt] options will be set automatically as well.",
294 "Again, if you think you know it better, please check the [TT]eigenfreq.xvg[tt]",
295 "output."
296 };
297
298 static gmx_bool bM = TRUE1, bCons = FALSE0;
299 static int begin = 1, end = 50, maxspec = 4000;
300 static real T = 298.15, width = 1;
301 t_pargs pa[] =
302 {
303 { "-m", FALSE0, etBOOL, {&bM},
304 "Divide elements of Hessian by product of sqrt(mass) of involved "
305 "atoms prior to diagonalization. This should be used for 'Normal Modes' "
306 "analysis" },
307 { "-first", FALSE0, etINT, {&begin},
308 "First eigenvector to write away" },
309 { "-last", FALSE0, etINT, {&end},
310 "Last eigenvector to write away" },
311 { "-maxspec", FALSE0, etINT, {&maxspec},
312 "Highest frequency (1/cm) to consider in the spectrum" },
313 { "-T", FALSE0, etREAL, {&T},
314 "Temperature for computing quantum heat capacity and enthalpy when using normal mode calculations to correct classical simulations" },
315 { "-constr", FALSE0, etBOOL, {&bCons},
316 "If constraints were used in the simulation but not in the normal mode analysis (this is the recommended way of doing it) you will need to set this for computing the quantum corrections." },
317 { "-width", FALSE0, etREAL, {&width},
318 "Width (sigma) of the gaussian peaks (1/cm) when generating a spectrum" }
319 };
320 FILE *out, *qc, *spec;
321 int status, trjout;
322 t_topology top;
323 gmx_mtop_t mtop;
324 int ePBC;
325 rvec *top_x;
326 matrix box;
327 real *eigenvalues;
328 real *eigenvectors;
329 real rdum, mass_fac, qcvtot, qutot, qcv, qu;
330 int natoms, ndim, nrow, ncol, count, nharm, nvsite;
331 char *grpname;
332 int i, j, k, l, d, gnx;
333 gmx_bool bSuck;
334 atom_id *index;
335 t_tpxheader tpx;
336 int version, generation;
337 real value, omega, nu;
338 real factor_gmx_to_omega2;
339 real factor_omega_to_wavenumber;
340 real *spectrum = NULL((void*)0);
341 real wfac;
342 output_env_t oenv;
343 const char *qcleg[] = {
344 "Heat Capacity cV (J/mol K)",
345 "Enthalpy H (kJ/mol)"
346 };
347 real * full_hessian = NULL((void*)0);
348 gmx_sparsematrix_t * sparse_hessian = NULL((void*)0);
349
350 t_filenm fnm[] = {
351 { efMTX, "-f", "hessian", ffREAD1<<1 },
352 { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
353 { efXVG, "-of", "eigenfreq", ffWRITE1<<2 },
354 { efXVG, "-ol", "eigenval", ffWRITE1<<2 },
355 { efXVG, "-os", "spectrum", ffOPTWR(1<<2| 1<<3) },
356 { efXVG, "-qc", "quant_corr", ffOPTWR(1<<2| 1<<3) },
357 { efTRN, "-v", "eigenvec", ffWRITE1<<2 }
358 };
359#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
360
361 if (!parse_common_args(&argc, argv, PCA_BE_NICE(1<<13),
1
Taking false branch
→
362 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
363 {
364 return 0;
365 }
366
367 /* Read tpr file for volume and number of harmonic terms */
368 read_tpxheader(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &tpx, TRUE1, &version, &generation);
369 snew(top_x, tpx.natoms)(top_x) = save_calloc("top_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmeig.c"
, 369, (tpx.natoms), sizeof(*(top_x)));
370
371 read_tpx(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0), box, &natoms,
372 top_x, NULL((void*)0), NULL((void*)0), &mtop);
373 if (bCons)
2
←
Taking false branch
→
374 {
375 nharm = get_nharm(&mtop, &nvsite);
376 }
377 else
378 {
379 nharm = 0;
380 nvsite = 0;
381 }
382 top = gmx_mtop_t_to_t_topology(&mtop);
383
384 bM = TRUE1;
385 ndim = DIM3*natoms;
386
387 if (opt2bSet("-qc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
3
←
Taking false branch
→
388 {
389 begin = 7+DIM3*nvsite;
390 end = DIM3*natoms;
391 }
392 if (begin < 1)
4
←
Taking false branch
→
393 {
394 begin = 1;
395 }
396 if (end > ndim)
5
←
Assuming 'end' is <= 'ndim'
→
6
←
Taking false branch
→
397 {
398 end = ndim;
399 }
400 printf("Using begin = %d and end = %d\n", begin, end);
401
402 /*open Hessian matrix */
403 gmx_mtxio_read(ftp2fn(efMTX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &nrow, &ncol, &full_hessian, &sparse_hessian);
7
←
Value assigned to 'sparse_hessian'
→
404
405 /* Memory for eigenvalues and eigenvectors (begin..end) */
406 snew(eigenvalues, nrow)(eigenvalues) = save_calloc("eigenvalues", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmeig.c"
, 406, (nrow), sizeof(*(eigenvalues)));
407 snew(eigenvectors, nrow*(end-begin+1))(eigenvectors) = save_calloc("eigenvectors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmeig.c"
, 407, (nrow*(end-begin+1)), sizeof(*(eigenvectors)));
408
409 /* If the Hessian is in sparse format we can calculate max (ndim-1) eigenvectors,
410 * and they must start at the first one. If this is not valid we convert to full matrix
411 * storage, but warn the user that we might run out of memory...
412 */
413 if ((sparse_hessian != NULL((void*)0)) && (begin != 1 || end == ndim))
8
←
Assuming 'sparse_hessian' is equal to null
→
414 {
415 if (begin != 1)
416 {
417 fprintf(stderrstderr, "Cannot use sparse Hessian with first eigenvector != 1.\n");
418 }
419 else if (end == ndim)
420 {
421 fprintf(stderrstderr, "Cannot use sparse Hessian to calculate all eigenvectors.\n");
422 }
423
424 fprintf(stderrstderr, "Will try to allocate memory and convert to full matrix representation...\n");
425
426 snew(full_hessian, nrow*ncol)(full_hessian) = save_calloc("full_hessian", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmeig.c"
, 426, (nrow*ncol), sizeof(*(full_hessian)));
427 for (i = 0; i < nrow*ncol; i++)
428 {
429 full_hessian[i] = 0;
430 }
431
432 for (i = 0; i < sparse_hessian->nrow; i++)
433 {
434 for (j = 0; j < sparse_hessian->ndata[i]; j++)
435 {
436 k = sparse_hessian->data[i][j].col;
437 value = sparse_hessian->data[i][j].value;
438 full_hessian[i*ndim+k] = value;
439 full_hessian[k*ndim+i] = value;
440 }
441 }
442 gmx_sparsematrix_destroy(sparse_hessian);
443 sparse_hessian = NULL((void*)0);
444 fprintf(stderrstderr, "Converted sparse to full matrix storage.\n");
445 }
446
447 if (full_hessian != NULL((void*)0))
9
←
Assuming 'full_hessian' is equal to null
→
10
←
Taking false branch
→
448 {
449 /* Using full matrix storage */
450 nma_full_hessian(full_hessian, nrow, bM, &top, begin, end,
451 eigenvalues, eigenvectors);
452 }
453 else
454 {
455 /* Sparse memory storage, allocate memory for eigenvectors */
456 snew(eigenvectors, ncol*end)(eigenvectors) = save_calloc("eigenvectors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmeig.c"
, 456, (ncol*end), sizeof(*(eigenvectors)));
457 nma_sparse_hessian(sparse_hessian, bM, &top, end, eigenvalues, eigenvectors);
11
←
Passing null pointer value via 1st parameter 'sparse_hessian'
→
12
←
Calling 'nma_sparse_hessian'
→
458 }
459
460 /* check the output, first 6 eigenvalues should be reasonably small */
461 bSuck = FALSE0;
462 for (i = begin-1; (i < 6); i++)
463 {
464 if (fabs(eigenvalues[i]) > 1.0e-3)
465 {
466 bSuck = TRUE1;
467 }
468 }
469 if (bSuck)
470 {
471 fprintf(stderrstderr, "\nOne of the lowest 6 eigenvalues has a non-zero value.\n");
472 fprintf(stderrstderr, "This could mean that the reference structure was not\n");
473 fprintf(stderrstderr, "properly energy minimized.\n");
474 }
475
476 /* now write the output */
477 fprintf (stderrstderr, "Writing eigenvalues...\n");
478 out = xvgropen(opt2fn("-ol", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
479 "Eigenvalues", "Eigenvalue index", "Eigenvalue [Gromacs units]",
480 oenv);
481 if (output_env_get_print_xvgr_codes(oenv))
482 {
483 if (bM)
484 {
485 fprintf(out, "@ subtitle \"mass weighted\"\n");
486 }
487 else
488 {
489 fprintf(out, "@ subtitle \"not mass weighted\"\n");
490 }
491 }
492
493 for (i = 0; i <= (end-begin); i++)
494 {
495 fprintf (out, "%6d %15g\n", begin+i, eigenvalues[i]);
496 }
497 gmx_ffclose(out);
498
499
500 if (opt2bSet("-qc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
501 {
502 qc = xvgropen(opt2fn("-qc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Quantum Corrections", "Eigenvector index", "", oenv);
503 xvgr_legend(qc, asize(qcleg)((int)(sizeof(qcleg)/sizeof((qcleg)[0]))), qcleg, oenv);
504 qcvtot = qutot = 0;
505 }
506 else
507 {
508 qc = NULL((void*)0);
509 }
510 printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
511
512 out = xvgropen(opt2fn("-of", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
513 "Eigenfrequencies", "Eigenvector index", "Wavenumber [cm\\S-1\\N]",
514 oenv);
515 if (output_env_get_print_xvgr_codes(oenv))
516 {
517 if (bM)
518 {
519 fprintf(out, "@ subtitle \"mass weighted\"\n");
520 }
521 else
522 {
523 fprintf(out, "@ subtitle \"not mass weighted\"\n");
524 }
525 }
526 /* Spectrum ? */
527 spec = NULL((void*)0);
528 if (opt2bSet("-os", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) && (maxspec > 0))
529 {
530 snew(spectrum, maxspec)(spectrum) = save_calloc("spectrum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_nmeig.c"
, 530, (maxspec), sizeof(*(spectrum)));
531 spec = xvgropen(opt2fn("-os", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
532 "Vibrational spectrum based on harmonic approximation",
533 "\\f{12}w\\f{4} (cm\\S-1\\N)",
534 "Intensity [Gromacs units]",
535 oenv);
536 for (i = 0; (i < maxspec); i++)
537 {
538 spectrum[i] = 0;
539 }
540 }
541
542 /* Gromacs units are kJ/(mol*nm*nm*amu),
543 * where amu is the atomic mass unit.
544 *
545 * For the eigenfrequencies we want to convert this to spectroscopic absorption
546 * wavenumbers given in cm^(-1), which is the frequency divided by the speed of
547 * light. Do this by first converting to omega^2 (units 1/s), take the square
548 * root, and finally divide by the speed of light (nm/ps in gromacs).
549 */
550 factor_gmx_to_omega2 = 1.0E21/(AVOGADRO(6.0221367e23)*AMU(1.6605402e-27));
551 factor_omega_to_wavenumber = 1.0E-5/(2.0*M_PI3.14159265358979323846*SPEED_OF_LIGHT(2.9979245800E05));
552
553 for (i = begin; (i <= end); i++)
554 {
555 value = eigenvalues[i-begin];
556 if (value < 0)
557 {
558 value = 0;
559 }
560 omega = sqrt(value*factor_gmx_to_omega2);
561 nu = 1e-12*omega/(2*M_PI3.14159265358979323846);
562 value = omega*factor_omega_to_wavenumber;
563 fprintf (out, "%6d %15g\n", i, value);
564 if (NULL((void*)0) != spec)
565 {
566 wfac = eigenvalues[i-begin]/(width*sqrt(2*M_PI3.14159265358979323846));
567 for (j = 0; (j < maxspec); j++)
568 {
569 spectrum[j] += wfac*exp(-sqr(j-value)/(2*sqr(width)));
570 }
571 }
572 if (NULL((void*)0) != qc)
573 {
574 qcv = cv_corr(nu, T);
575 qu = u_corr(nu, T);
576 if (i > end-nharm)
577 {
578 qcv += BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*KILO(1e3);
579 qu += BOLTZ(((1.380658e-23)*(6.0221367e23))/(1e3))*T;
580 }
581 fprintf (qc, "%6d %15g %15g\n", i, qcv, qu);
582 qcvtot += qcv;
583 qutot += qu;
584 }
585 }
586 gmx_ffclose(out);
587 if (NULL((void*)0) != spec)
588 {
589 for (j = 0; (j < maxspec); j++)
590 {
591 fprintf(spec, "%10g %10g\n", 1.0*j, spectrum[j]);
592 }
593 gmx_ffclose(spec);
594 }
595 if (NULL((void*)0) != qc)
596 {
597 printf("Quantum corrections for harmonic degrees of freedom\n");
598 printf("Use appropriate -first and -last options to get reliable results.\n");
599 printf("There were %d constraints and %d vsites in the simulation\n",
600 nharm, nvsite);
601 printf("Total correction to cV = %g J/mol K\n", qcvtot);
602 printf("Total correction to H = %g kJ/mol\n", qutot);
603 gmx_ffclose(qc);
604 please_cite(stdoutstdout, "Caleman2011b");
605 }
606 /* Writing eigenvectors. Note that if mass scaling was used, the eigenvectors
607 * were scaled back from mass weighted cartesian to plain cartesian in the
608 * nma_full_hessian() or nma_sparse_hessian() routines. Mass scaled vectors
609 * will not be strictly orthogonal in plain cartesian scalar products.
610 */
611 write_eigenvectors(opt2fn("-v", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), natoms, eigenvectors, FALSE0, begin, end,
612 eWXR_NO, NULL((void*)0), FALSE0, top_x, bM, eigenvalues);
613
614 return 0;
615}