src/gromacs/gmxana/gmx_tcaf.cpp

   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,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source 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 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstdio>
  42 #include <cstring>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/commandline/viewit.h"
  46 #include "gromacs/correlationfunctions/autocorr.h"
  47 #include "gromacs/correlationfunctions/expfit.h"
  48 #include "gromacs/fileio/confio.h"
  49 #include "gromacs/fileio/trxio.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/gmxana/gstat.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/topology/index.h"
  58 #include "gromacs/topology/topology.h"
  59 #include "gromacs/trajectory/trajectoryframe.h"
  60 #include "gromacs/utility/arraysize.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/futil.h"
  63 #include "gromacs/utility/gmxassert.h"
  64 #include "gromacs/utility/smalloc.h"
  65
  66 #define NK 24
  67 #define NPK 4
  68
  69 #define NKC 6
  70 #define NKC0 4
  71 static const int kset_c[NKC + 1] = { 0, 3, 9, 13, 16, 19, NK };
  72
  73 static rvec v0[NK] = { { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },  { 1, 1, 0 },  { 1, -1, 0 },
  74                        { 1, 0, 1 },  { 1, 0, -1 }, { 0, 1, 1 },  { 0, 1, -1 }, { 1, 1, 1 },
  75                        { 1, 1, -1 }, { 1, -1, 1 }, { -1, 1, 1 }, { 2, 0, 0 },  { 0, 2, 0 },
  76                        { 0, 0, 2 },  { 3, 0, 0 },  { 0, 3, 0 },  { 0, 0, 3 },  { 4, 0, 0 },
  77                        { 0, 4, 0 },  { 0, 0, 4 } };
  78 static rvec v1[NK] = { { 0, 1, 0 },  { 0, 0, 1 },  { 1, 0, 0 },  { 0, 0, 1 }, { 0, 0, 1 },
  79                        { 0, 1, 0 },  { 0, 1, 0 },  { 1, 0, 0 },  { 1, 0, 0 }, { 1, -1, 0 },
  80                        { 1, -1, 0 }, { 1, 0, -1 }, { 0, 1, -1 }, { 0, 1, 0 }, { 0, 0, 1 },
  81                        { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },  { 1, 0, 0 }, { 0, 1, 0 },
  82                        { 0, 0, 1 },  { 1, 0, 0 } };
  83 static rvec v2[NK] = { { 0, 0, 1 },  { 1, 0, 0 }, { 0, 1, 0 },  { 1, -1, 0 }, { 1, 1, 0 },
  84                        { 1, 0, -1 }, { 1, 0, 1 }, { 0, 1, -1 }, { 0, 1, 1 },  { 1, 1, -2 },
  85                        { 1, 1, 2 },  { 1, 2, 1 }, { 2, 1, 1 },  { 0, 0, 1 },  { 1, 0, 0 },
  86                        { 0, 1, 0 },  { 0, 0, 1 }, { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },
  87                        { 1, 0, 0 },  { 0, 1, 0 } };
  88
  89 static void process_tcaf(int                     nframes,
  90                          real                    dt,
  91                          int                     nkc,
  92                          real**                  tc,
  93                          rvec*                   kfac,
  94                          real                    rho,
  95                          real                    wt,
  96                          const char*             fn_trans,
  97                          const char*             fn_tca,
  98                          const char*             fn_tc,
  99                          const char*             fn_tcf,
 100                          const char*             fn_cub,
 101                          const char*             fn_vk,
 102                          const gmx_output_env_t* oenv)
 103 {
 104     FILE * fp, *fp_vk, *fp_cub = nullptr;
 105     int    nk, ntc;
 106     real **tcaf, **tcafc = nullptr, eta, *sig;
 107     int    i, j, k, kc;
 108     int    ncorr;
 109     double fitparms[3];
 110
 111     nk  = kset_c[nkc];
 112     ntc = nk * NPK;
 113
 114     if (fn_trans)
 115     {
 116         fp = xvgropen(fn_trans, "Transverse Current", "Time (ps)", "TC (nm/ps)", oenv);
 117         for (i = 0; i < nframes; i++)
 118         {
 119             fprintf(fp, "%g", i * dt);
 120             for (j = 0; j < ntc; j++)
 121             {
 122                 fprintf(fp, " %g", tc[j][i]);
 123             }
 124             fprintf(fp, "\n");
 125         }
 126         xvgrclose(fp);
 127         do_view(oenv, fn_trans, "-nxy");
 128     }
 129
 130     ncorr = (nframes + 1) / 2;
 131     if (ncorr > gmx::roundToInt(5 * wt / dt))
 132     {
 133         ncorr = gmx::roundToInt(5 * wt / dt) + 1;
 134     }
 135     snew(tcaf, nk);
 136     for (k = 0; k < nk; k++)
 137     {
 138         snew(tcaf[k], ncorr);
 139     }
 140     if (fn_cub)
 141     {
 142         snew(tcafc, nkc);
 143         for (k = 0; k < nkc; k++)
 144         {
 145             snew(tcafc[k], ncorr);
 146         }
 147     }
 148     snew(sig, ncorr);
 149     for (i = 0; i < ncorr; i++)
 150     {
 151         sig[i] = std::exp(0.5 * i * dt / wt);
 152     }
 153
 154     low_do_autocorr(fn_tca,
 155                     oenv,
 156                     "Transverse Current Autocorrelation Functions",
 157                     nframes,
 158                     ntc,
 159                     ncorr,
 160                     tc,
 161                     dt,
 162                     eacNormal,
 163                     1,
 164                     FALSE,
 165                     FALSE,
 166                     FALSE,
 167                     0,
 168                     0,
 169                     0);
 170     do_view(oenv, fn_tca, "-nxy");
 171
 172     fp = xvgropen(fn_tc, "Transverse Current Autocorrelation Functions", "Time (ps)", "TCAF", oenv);
 173     for (i = 0; i < ncorr; i++)
 174     {
 175         kc = 0;
 176         fprintf(fp, "%g", i * dt);
 177         for (k = 0; k < nk; k++)
 178         {
 179             for (j = 0; j < NPK; j++)
 180             {
 181                 tcaf[k][i] += tc[NPK * k + j][i];
 182             }
 183             if (fn_cub)
 184             {
 185                 for (j = 0; j < NPK; j++)
 186                 {
 187                     tcafc[kc][i] += tc[NPK * k + j][i];
 188                 }
 189             }
 190             if (i == 0)
 191             {
 192                 fprintf(fp, " %g", 1.0);
 193             }
 194             else
 195             {
 196                 tcaf[k][i] /= tcaf[k][0];
 197                 fprintf(fp, " %g", tcaf[k][i]);
 198             }
 199             if (k + 1 == kset_c[kc + 1])
 200             {
 201                 kc++;
 202             }
 203         }
 204         fprintf(fp, "\n");
 205     }
 206     xvgrclose(fp);
 207     do_view(oenv, fn_tc, "-nxy");
 208
 209     if (fn_cub)
 210     {
 211         fp_cub = xvgropen(fn_cub, "TCAFs and fits", "Time (ps)", "TCAF", oenv);
 212         for (kc = 0; kc < nkc; kc++)
 213         {
 214             fprintf(fp_cub, "%g %g\n", 0.0, 1.0);
 215             for (i = 1; i < ncorr; i++)
 216             {
 217                 tcafc[kc][i] /= tcafc[kc][0];
 218                 fprintf(fp_cub, "%g %g\n", i * dt, tcafc[kc][i]);
 219             }
 220             fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 221             tcafc[kc][0] = 1.0;
 222         }
 223     }
 224
 225     fp_vk = xvgropen(fn_vk, "Fits", "k (nm\\S-1\\N)", "\\8h\\4 (10\\S-3\\N kg m\\S-1\\N s\\S-1\\N)", oenv);
 226     if (output_env_get_print_xvgr_codes(oenv))
 227     {
 228         fprintf(fp_vk, "@    s0 symbol 2\n");
 229         fprintf(fp_vk, "@    s0 symbol color 1\n");
 230         fprintf(fp_vk, "@    s0 linestyle 0\n");
 231         if (fn_cub)
 232         {
 233             fprintf(fp_vk, "@    s1 symbol 3\n");
 234             fprintf(fp_vk, "@    s1 symbol color 2\n");
 235         }
 236     }
 237     fp = xvgropen(fn_tcf, "TCAF Fits", "Time (ps)", "", oenv);
 238     for (k = 0; k < nk; k++)
 239     {
 240         tcaf[k][0]  = 1.0;
 241         fitparms[0] = 1;
 242         fitparms[1] = 1;
 243         do_lmfit(ncorr, tcaf[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr);
 244         eta = 1000 * fitparms[1] * rho
 245               / (4 * fitparms[0] * gmx::c_pico * norm2(kfac[k]) / (gmx::c_nano * gmx::c_nano));
 246         fprintf(stdout, "k %6.3f  tau %6.3f  eta %8.5f 10^-3 kg/(m s)\n", norm(kfac[k]), fitparms[0], eta);
 247         fprintf(fp_vk, "%6.3f %g\n", norm(kfac[k]), eta);
 248         for (i = 0; i < ncorr; i++)
 249         {
 250             fprintf(fp, "%g %g\n", i * dt, fit_function(effnVAC, fitparms, i * dt));
 251         }
 252         fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 253     }
 254     xvgrclose(fp);
 255     do_view(oenv, fn_tcf, "-nxy");
 256
 257     if (fn_cub)
 258     {
 259         fprintf(stdout, "Averaged over k-vectors:\n");
 260         fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 261         for (k = 0; k < nkc; k++)
 262         {
 263             tcafc[k][0] = 1.0;
 264             fitparms[0] = 1;
 265             fitparms[1] = 1;
 266             do_lmfit(ncorr, tcafc[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr);
 267             eta = 1000 * fitparms[1] * rho
 268                   / (4 * fitparms[0] * gmx::c_pico * norm2(kfac[kset_c[k]]) / (gmx::c_nano * gmx::c_nano));
 269             fprintf(stdout,
 270                     "k %6.3f  tau %6.3f  Omega %6.3f  eta %8.5f 10^-3 kg/(m s)\n",
 271                     norm(kfac[kset_c[k]]),
 272                     fitparms[0],
 273                     fitparms[1],
 274                     eta);
 275             fprintf(fp_vk, "%6.3f %g\n", norm(kfac[kset_c[k]]), eta);
 276             for (i = 0; i < ncorr; i++)
 277             {
 278                 fprintf(fp_cub, "%g %g\n", i * dt, fit_function(effnVAC, fitparms, i * dt));
 279             }
 280             fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 281         }
 282         fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 283         xvgrclose(fp_cub);
 284         do_view(oenv, fn_cub, "-nxy");
 285     }
 286     xvgrclose(fp_vk);
 287     do_view(oenv, fn_vk, "-nxy");
 288 }
 289
 290
 291 int gmx_tcaf(int argc, char* argv[])
 292 {
 293     const char* desc[] = {
 294         "[THISMODULE] computes tranverse current autocorrelations.",
 295         "These are used to estimate the shear viscosity, [GRK]eta[grk].",
 296         "For details see: Palmer, Phys. Rev. E 49 (1994) pp 359-366.[PAR]",
 297         "Transverse currents are calculated using the",
 298         "k-vectors (1,0,0) and (2,0,0) each also in the [IT]y[it]- and [IT]z[it]-direction,",
 299         "(1,1,0) and (1,-1,0) each also in the 2 other planes (these vectors",
 300         "are not independent) and (1,1,1) and the 3 other box diagonals (also",
 301         "not independent). For each k-vector the sine and cosine are used, in",
 302         "combination with the velocity in 2 perpendicular directions. This gives",
 303         "a total of 16*2*2=64 transverse currents. One autocorrelation is",
 304         "calculated fitted for each k-vector, which gives 16 TCAFs. Each of",
 305         "these TCAFs is fitted to [MATH]f(t) = [EXP]-v[exp]([COSH]Wv[cosh] + 1/W ",
 306         "[SINH]Wv[sinh])[math],",
 307         "[MATH]v = -t/(2 [GRK]tau[grk])[math], [MATH]W = [SQRT]1 - 4 [GRK]tau[grk] ",
 308         "[GRK]eta[grk]/[GRK]rho[grk] k^2[sqrt][math], which gives 16 values of [GRK]tau[grk]",
 309         "and [GRK]eta[grk]. The fit weights decay exponentially with time constant [MATH]w[math] ",
 310         "(given with [TT]-wt[tt]) as [MATH][EXP]-t/w[exp][math], and the TCAF and",
 311         "fit are calculated up to time [MATH]5*w[math].",
 312         "The [GRK]eta[grk] values should be fitted to [MATH]1 - a [GRK]eta[grk](k) k^2[math], ",
 313         "from which one can estimate the shear viscosity at k=0.[PAR]",
 314         "When the box is cubic, one can use the option [TT]-oc[tt], which",
 315         "averages the TCAFs over all k-vectors with the same length.",
 316         "This results in more accurate TCAFs.",
 317         "Both the cubic TCAFs and fits are written to [TT]-oc[tt]",
 318         "The cubic [GRK]eta[grk] estimates are also written to [TT]-ov[tt].[PAR]",
 319         "With option [TT]-mol[tt], the transverse current is determined of",
 320         "molecules instead of atoms. In this case, the index group should",
 321         "consist of molecule numbers instead of atom numbers.[PAR]",
 322         "The k-dependent viscosities in the [TT]-ov[tt] file should be",
 323         "fitted to [MATH][GRK]eta[grk](k) = [GRK]eta[grk][SUB]0[sub] (1 - a k^2)[math] to obtain ",
 324         "the viscosity at",
 325         "infinite wavelength.[PAR]",
 326         "[BB]Note:[bb] make sure you write coordinates and velocities often enough.",
 327         "The initial, non-exponential, part of the autocorrelation function",
 328         "is very important for obtaining a good fit."
 329     };
 330
 331     static gmx_bool bMol = FALSE, bK34 = FALSE;
 332     static real     wt   = 5;
 333     t_pargs         pa[] = {
 334         { "-mol", FALSE, etBOOL, { &bMol }, "Calculate TCAF of molecules" },
 335         { "-k34", FALSE, etBOOL, { &bK34 }, "Also use k=(3,0,0) and k=(4,0,0)" },
 336         { "-wt", FALSE, etREAL, { &wt }, "Exponential decay time for the TCAF fit weights" }
 337     };
 338
 339     t_topology        top;
 340     PbcType           pbcType;
 341     t_trxframe        fr;
 342     matrix            box;
 343     gmx_bool          bTop;
 344     int               gnx;
 345     int *             index, *atndx = nullptr, at;
 346     char*             grpname;
 347     char              title[256];
 348     real              t0, t1, dt, m, mtot, sysmass, rho, sx, cx;
 349     t_trxstatus*      status;
 350     int               nframes, n_alloc, i, j, k, d;
 351     rvec              mv_mol, cm_mol, kfac[NK];
 352     int               nkc, nk, ntc;
 353     real**            tc;
 354     gmx_output_env_t* oenv;
 355
 356 #define NHISTO 360
 357
 358     t_filenm fnm[] = { { efTRN, "-f", nullptr, ffREAD },      { efTPS, nullptr, nullptr, ffOPTRD },
 359                        { efNDX, nullptr, nullptr, ffOPTRD },  { efXVG, "-ot", "transcur", ffOPTWR },
 360                        { efXVG, "-oa", "tcaf_all", ffWRITE }, { efXVG, "-o", "tcaf", ffWRITE },
 361                        { efXVG, "-of", "tcaf_fit", ffWRITE }, { efXVG, "-oc", "tcaf_cub", ffOPTWR },
 362                        { efXVG, "-ov", "visc_k", ffWRITE } };
 363 #define NFILE asize(fnm)
 364     int      npargs;
 365     t_pargs* ppa;
 366
 367     npargs = asize(pa);
 368     ppa    = add_acf_pargs(&npargs, pa);
 369
 370     if (!parse_common_args(
 371                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
 372     {
 373         sfree(ppa);
 374         return 0;
 375     }
 376
 377     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
 378     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 379
 380     if (bMol)
 381     {
 382         if (!bTop)
 383         {
 384             gmx_fatal(FARGS, "Need a topology to determine the molecules");
 385         }
 386         atndx = top.mols.index;
 387     }
 388
 389     if (bK34)
 390     {
 391         nkc = NKC;
 392     }
 393     else
 394     {
 395         nkc = NKC0;
 396     }
 397     nk = kset_c[nkc];
 398     GMX_ASSERT(nk >= 16, "Has to be over 16 because nkc is either NKC or NKC0.");
 399     ntc = nk * NPK;
 400
 401     sprintf(title, "Velocity Autocorrelation Function for %s", grpname);
 402
 403     sysmass = 0;
 404     for (i = 0; i < nk; i++)
 405     {
 406         if (iprod(v0[i], v1[i]) != 0)
 407         {
 408             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
 409         }
 410         if (iprod(v0[i], v2[i]) != 0)
 411         {
 412             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
 413         }
 414         if (iprod(v1[i], v2[i]) != 0)
 415         {
 416             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
 417         }
 418         unitv(v1[i], v1[i]);
 419         unitv(v2[i], v2[i]);
 420     }
 421     snew(tc, ntc);
 422     for (i = 0; i < top.atoms.nr; i++)
 423     {
 424         sysmass += top.atoms.atom[i].m;
 425     }
 426
 427     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_X | TRX_NEED_V);
 428     t0 = fr.time;
 429
 430     n_alloc = 0;
 431     nframes = 0;
 432     rho     = 0;
 433
 434     do
 435     {
 436
 437         if (nframes >= n_alloc)
 438         {
 439             n_alloc += 100;
 440             for (i = 0; i < ntc; i++)
 441             {
 442                 srenew(tc[i], n_alloc);
 443             }
 444         }
 445
 446         rho += 1 / det(fr.box);
 447         for (k = 0; k < nk; k++)
 448         {
 449             for (d = 0; d < DIM; d++)
 450             {
 451                 kfac[k][d] = 2 * M_PI * v0[k][d] / fr.box[d][d];
 452             }
 453         }
 454         for (i = 0; i < ntc; i++)
 455         {
 456             tc[i][nframes] = 0;
 457         }
 458
 459         for (i = 0; i < gnx; i++)
 460         {
 461             if (bMol)
 462             {
 463                 clear_rvec(mv_mol);
 464                 clear_rvec(cm_mol);
 465                 mtot = 0;
 466                 for (j = 0; j < atndx[index[i] + 1] - atndx[index[i]]; j++)
 467                 {
 468                     at = atndx[index[i]] + j;
 469                     m  = top.atoms.atom[at].m;
 470                     mv_mol[XX] += m * fr.v[at][XX];
 471                     mv_mol[YY] += m * fr.v[at][YY];
 472                     mv_mol[ZZ] += m * fr.v[at][ZZ];
 473                     cm_mol[XX] += m * fr.x[at][XX];
 474                     cm_mol[YY] += m * fr.x[at][YY];
 475                     cm_mol[ZZ] += m * fr.x[at][ZZ];
 476                     mtot += m;
 477                 }
 478                 svmul(1.0 / mtot, cm_mol, cm_mol);
 479             }
 480             else
 481             {
 482                 svmul(top.atoms.atom[index[i]].m, fr.v[index[i]], mv_mol);
 483             }
 484
 485             if (!bMol)
 486             {
 487                 copy_rvec(fr.x[index[i]], cm_mol);
 488             }
 489             j = 0;
 490             for (k = 0; k < nk; k++)
 491             {
 492                 sx = std::sin(iprod(kfac[k], cm_mol));
 493                 cx = std::cos(iprod(kfac[k], cm_mol));
 494                 tc[j][nframes] += sx * iprod(v1[k], mv_mol);
 495                 j++;
 496                 tc[j][nframes] += cx * iprod(v1[k], mv_mol);
 497                 j++;
 498                 tc[j][nframes] += sx * iprod(v2[k], mv_mol);
 499                 j++;
 500                 tc[j][nframes] += cx * iprod(v2[k], mv_mol);
 501                 j++;
 502             }
 503         }
 504
 505         t1 = fr.time;
 506         nframes++;
 507     } while (read_next_frame(oenv, status, &fr));
 508     close_trx(status);
 509
 510     dt = (t1 - t0) / (nframes - 1);
 511
 512     rho *= sysmass / nframes * gmx::c_amu / (gmx::c_nano * gmx::c_nano * gmx::c_nano);
 513     fprintf(stdout, "Density = %g (kg/m^3)\n", rho);
 514     process_tcaf(nframes,
 515                  dt,
 516                  nkc,
 517                  tc,
 518                  kfac,
 519                  rho,
 520                  wt,
 521                  opt2fn_null("-ot", NFILE, fnm),
 522                  opt2fn("-oa", NFILE, fnm),
 523                  opt2fn("-o", NFILE, fnm),
 524                  opt2fn("-of", NFILE, fnm),
 525                  opt2fn_null("-oc", NFILE, fnm),
 526                  opt2fn("-ov", NFILE, fnm),
 527                  oenv);
 528
 529     return 0;
 530 }