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