src/tools/gmx_dos.c

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