C++ opt2fn函数代码示例

int gmx_gmx2pqr(int argc, char *argv[])
{
    const char *desc[] = {
        "\tThis convert from gromacs to APBS pqr files.",
        "Use the -a1 flag for the start of the bond vector",
        "and the -a2 flag for the end of the bond vector.",
        "\n\tUse the -exclude flag to specify atoms which you",
        "would like to know the speficic contribution to the",
        "electrostatic field at the bond midpoint, or would",
        "like to exclude from the field calculation.",
        "Use the -site flag to add dummy atoms like those used",
        "by the Cho group {Choi, J. H.; Oh, K. I.; Lee, H.; Lee,",
        "C.; Cho, M.; Nitrile and thiocyanate IR",
        "probes: Quantum chemistry calculation studies and",
        "multivariate least-square fitting analysis; J. Phys.",
        "Chem. 2008}.  By default, both -site and -exclude flag",
        "will include the atoms indicated by -a1 and -a2.",
        "\n\tUse the -select flag to specify which atoms you",
        "want to include in the .pqr file and all calculations.",
        "For example, '(not resname SOL and not name Na) or",
        "resname SOL and same residue as within 0.5 of resname",
        "CNC and (name NE or name CD or name CB or name SG)",
        "would give protein and water near the heavy backbone",
        "atoms of any cyanocysteine residues. By default, when",
        "the -exclude flag is used, potentials will be reported",
        "in two group--water (resname HOH or resname SOL) and",
        "non-water contributions."
        "\n"
    };
    /* Command-line arguments */
    int             a1 = 0; // initialized for error checking
    int             a2 = 0;
    real            pdie = 1;
    real            delta = 0.010; // nm
    real            ring_dist = 0.07; // nm
    static char     *exclude = NULL;
    static char     *site = NULL;
    bool            bVerbose = false; // 0=false, 1=true
    bool            dopqr = true;
    
    t_pargs         pa[] = {
        { "-a1", TRUE, etINT,
            {&a1}, "Starting atom for bond vector--ie: CD in CNC"},
        { "-a2", TRUE, etINT,
            {&a2}, "Ending atom for bond vector--ie: NE in CNC"},
        { "-site", TRUE, etSTR,
            {&site}, "Atoms to be included as sites for the electrostatic model used by the Cho group.  Atoms pass by -a1 and -a2 are added to this group automatically.  Note that the whole selection string will need to be quoted so that your shell will pass it in as a string."},
        { "-exclude", TRUE, etSTR,
            {&exclude}, "Specific atoms to exclude from the field calculations.  Atoms pass by -a1 and -a2 are added to this group automatically.  Note that the whole selection string will need to be quoted so that your shell will pass it in as a string."},
        { "-pdie", TRUE, etREAL,
            {&pdie}, "Protein dielectric constant"},
        { "-delta", TRUE, etREAL,
            {&delta}, "(nm) Spacing between dummy atoms at the bond midpoint"},
        { "-dopqr", TRUE, etBOOL,
            {&dopqr}, "Write out a .pqr file.  Default: True"},
        { "-v", FALSE, etBOOL,
            {&bVerbose}, "Be slightly more verbose"}
    };
    
    t_filenm        fnm[] = {
         { efDAT, "-d", "/Users/ritchie/Utilities/apbs/AMBER.DAT", ffREAD }
        ,{ efPQR, "-o", NULL, ffWRITE }
        ,{ efXVG, "-of", "coulomb_field", ffWRITE }
        ,{ efDAT, "-bw", "bw.dat", ffOPTWR}
        ,{ efXVG, "-or", "coulomb_field_per_residue", ffOPTWR}
        //,{ efNDX, NULL, NULL, ffREAD }
    };
    
#define NFILE asize(fnm)
    
    gmx_ana_traj_t      *trj;
    t_topology          *top;
    output_env_t        oenv;
    //t_dsdata            d;
    //t_analysisdata d;
    t_analysisdata      d;
    int                 ngrps;
    gmx_ana_selection_t **sel;
    int                 g;
    int                 rc;
    
    CopyRight(stderr, argv[0]);
    
    gmx_ana_traj_create(&trj, ANA_REQUIRE_TOP);
    gmx_ana_set_nanagrps(trj, -1);
    parse_trjana_args(trj, &argc, argv, 0,
                      NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
                      &oenv);
    gmx_ana_get_topology(trj, FALSE, &top, NULL);

    gmx_ana_get_nanagrps(trj, &ngrps);
    gmx_ana_get_anagrps(trj, &sel);
    gmx_ana_init_coverfrac(trj, CFRAC_SOLIDANGLE);
    
    /* Get output file names */
    const char          *datname, *framepqr;
    datname = opt2fn("-d", NFILE, fnm);
    framepqr = opt2fn("-o", NFILE, fnm);
    
    /* open xvg file */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    fprintf(stderr, "\nWARNING: The simulation uses %s temperature and/or %s pressure coupling,\n"
                            "         the continuation will only be exact when an energy file is supplied\n\n",
                            ETCOUPLTYPE(etcNOSEHOOVER),
                            EPCOUPLTYPE(epcPARRINELLORAHMAN));
                }
            }
            if (bFepState)
            {
                ir->fepvals->init_fep_state = init_fep_state;
            }
        }
    }

    if (bNsteps)
    {
        fprintf(stderr, "Setting nsteps to %s\n", gmx_step_str(nsteps_req, buf));
        ir->nsteps = nsteps_req;
    }
    else
    {
        /* Determine total number of steps remaining */
        if (bExtend)
        {
            ir->nsteps = ir->nsteps - (run_step - ir->init_step) + (gmx_int64_t)(extend_t/ir->delta_t + 0.5);
            printf("Extending remaining runtime of by %g ps (now %s steps)\n",
                   extend_t, gmx_step_str(ir->nsteps, buf));
        }
        else if (bUntil)
        {
            printf("nsteps = %s, run_step = %s, current_t = %g, until = %g\n",
                   gmx_step_str(ir->nsteps, buf),
                   gmx_step_str(run_step, buf2),
                   run_t, until_t);
            ir->nsteps = (gmx_int64_t)((until_t - run_t)/ir->delta_t + 0.5);
            printf("Extending remaining runtime until %g ps (now %s steps)\n",
                   until_t, gmx_step_str(ir->nsteps, buf));
        }
        else
        {
            ir->nsteps -= run_step - ir->init_step;
            /* Print message */
            printf("%s steps (%g ps) remaining from first run.\n",
                   gmx_step_str(ir->nsteps, buf), ir->nsteps*ir->delta_t);
        }
    }

    if (bNsteps || bZeroQ || (ir->nsteps > 0))
    {
        ir->init_step = run_step;

        if (ftp2bSet(efNDX, NFILE, fnm) ||
            !(bNsteps || bExtend || bUntil || bTraj))
        {
            atoms = gmx_mtop_global_atoms(&mtop);
            get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1,
                      &gnx, &index, &grpname);
            if (!bZeroQ)
            {
                bSel = (gnx != state.natoms);
                for (i = 0; ((i < gnx) && (!bSel)); i++)
                {
                    bSel = (i != index[i]);
                }
            }
            else
            {
                bSel = FALSE;
            }
            if (bSel)
            {
                fprintf(stderr, "Will write subset %s of original tpx containing %d "
                        "atoms\n", grpname, gnx);
                reduce_topology_x(gnx, index, &mtop, state.x, state.v);
                state.natoms = gnx;
            }
            else if (bZeroQ)
            {
                zeroq(index, &mtop);
                fprintf(stderr, "Zero-ing charges for group %s\n", grpname);
            }
            else
            {
                fprintf(stderr, "Will write full tpx file (no selection)\n");
            }
        }

        state_t = ir->init_t + ir->init_step*ir->delta_t;
        sprintf(buf,   "Writing statusfile with starting step %s%s and length %s%s steps...\n", "%10", GMX_PRId64, "%10", GMX_PRId64);
        fprintf(stderr, buf, ir->init_step, ir->nsteps);
        fprintf(stderr, "                                 time %10.3f and length %10.3f ps\n",
                state_t, ir->nsteps*ir->delta_t);
        write_tpx_state(opt2fn("-o", NFILE, fnm), ir, &state, &mtop);
    }
    else
    {
        printf("You've simulated long enough. Not writing tpr file\n");
    }

    return 0;
}

//.........这里部分代码省略.........
	    npcheck++;
	}
      }
    }
    if (npcheck != nphobic)
      gmx_fatal(FARGS,"Consistency check failed: not all %d atoms in the hydrophobic index\n"
		  "found in the normal index selection (%d atoms)",nphobic,npcheck);
  }
  else
    nphobic = 0;
    
  for(i=0; (i<nx[0]); i++) {
    ii = index[0][i];
    if (!bFindex) {
      bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;
      if (bPhobic[i] && bOut[ii])
	nphobic++;
    }
    if (bDGsol)
      if (!gmx_atomprop_query(aps,epropDGsol,
			      *(atoms->resinfo[atoms->atom[ii].resind].name),
			      *(atoms->atomtype[ii]),&(dgs_factor[i])))
	dgs_factor[i] = dgs_default;
    if (debug)
      fprintf(debug,"Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s\n",
	      ii+1,*(atoms->resinfo[atoms->atom[ii].resind].name),
	      *(atoms->atomname[ii]),
	      atoms->atom[ii].q,radius[ii]-solsize,dgs_factor[i],
	      BOOL(bPhobic[i]));
  }
  fprintf(stderr,"%d out of %d atoms were classified as hydrophobic\n",
	  nphobic,nx[1]);
  
  fp=xvgropen(opt2fn("-o",nfile,fnm),"Solvent Accessible Surface","Time (ps)",
	      "Area (nm\\S2\\N)",oenv);
  xvgr_legend(fp,asize(flegend) - (bDGsol ? 0 : 1),flegend,oenv);
  vfile = opt2fn_null("-tv",nfile,fnm);
  if (vfile) {
    if (!bTop) {
      gmx_fatal(FARGS,"Need a tpr file for option -tv");
    }
    vp=xvgropen(vfile,"Volume and Density","Time (ps)","",oenv);
    xvgr_legend(vp,asize(vlegend),vlegend,oenv);
    totmass  = 0;
    ndefault = 0;
    for(i=0; (i<nx[0]); i++) {
      real mm;
      ii = index[0][i];
      /*
      if (!query_atomprop(atomprop,epropMass,
			  *(top->atoms.resname[top->atoms.atom[ii].resnr]),
			  *(top->atoms.atomname[ii]),&mm))
	ndefault++;
      totmass += mm;
      */
      totmass += atoms->atom[ii].m;
    }
    if (ndefault)
      fprintf(stderr,"WARNING: Using %d default masses for density calculation, which most likely are inaccurate\n",ndefault);
  }
  else
    vp = NULL;
    
  gmx_atomprop_destroy(aps);

  if (bPBC)

//.........这里部分代码省略.........
    const char  *TargetFile;
    const char  *OriginFile;
    const char  *EigvecFile;

    output_env_t oenv;

    /*to read topology file*/
    t_topology  top;
    int         ePBC;
    char        title[STRLEN];
    matrix      topbox;
    rvec       *xtop;
    gmx_bool    bTop, bFit1;

    t_filenm    fnm[] = {
        { efTRN, "-f",    "eigenvec",    ffREAD  },
        { efXVG, "-eig",  "eigenval",    ffOPTRD },
        { efTPS, NULL,    NULL,          ffREAD },
        { efNDX, NULL,    NULL,  ffOPTRD },
        { efSTX, "-tar", "target", ffOPTRD},
        { efSTX, "-ori", "origin", ffOPTRD},
        { efEDI, "-o", "sam", ffWRITE }
    };
#define NFILE asize(fnm)
    edi_params.outfrq = 100; edi_params.slope = 0.0; edi_params.maxedsteps = 0;
    if (!parse_common_args(&argc, argv, 0,
                           NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }

    indexfile       = ftp2fn_null(efNDX, NFILE, fnm);
    EdiFile         = ftp2fn(efEDI, NFILE, fnm);
    TargetFile      = opt2fn_null("-tar", NFILE, fnm);
    OriginFile      = opt2fn_null("-ori", NFILE, fnm);


    for (ev_class = 0; ev_class < evNr; ++ev_class)
    {
        if (opt2parg_bSet(evOptions[ev_class], NPA, pa))
        {
            /*get list of eigenvectors*/
            nvecs = sscan_list(&(listen[ev_class]), opt2parg_str(evOptions[ev_class], NPA, pa), evOptions[ev_class]);
            if (ev_class < evStepNr-2)
            {
                /*if apropriate get list of stepsizes for these eigenvectors*/
                if (opt2parg_bSet(evStepOptions[ev_class], NPA, pa))
                {
                    evStepList[ev_class] =
                        scan_vecparams(opt2parg_str(evStepOptions[ev_class], NPA, pa), evStepOptions[ev_class], nvecs);
                }
                else   /*if list is not given fill with zeros */
                {
                    snew(evStepList[ev_class], nvecs);
                    for (i = 0; i < nvecs; i++)
                    {
                        evStepList[ev_class][i] = 0.0;
                    }
                }
            }
            else if (ev_class == evRADFIX)
            {
                snew(evStepList[ev_class], nvecs);
                for (i = 0; i < nvecs; i++)
                {
                    evStepList[ev_class][i] = radstep;

int main(int argc,char *argv[])
{
  static char *desc[] = {
    "compnl compares two neighborlists as generated by [TT]mdrun[tt]",
    "in the log file, when the environment variable DUMPNL is set to",
    "a number larger than 0. compnl is mainly used for debugging the",
    "mdrun internals and not for end-users."
  };
  FILE    *in,*out;
  int     i,j,nmiss,mod;
  char    **fn,title[256];
  int     ***mat,nnb;
  real    mb;
  bool    bConf;
  rvec    *x = NULL;
  rvec    dx;
  matrix  box;
  t_atoms atoms;
  t_pbc   pbc;
  
  t_filenm fnm[] = {
    { efLOG, "-f1", NULL, ffREAD },
    { efLOG, "-f2", NULL, ffREAD },
    { efOUT, "-o",  "compnl", ffWRITE },
    { efSTX, "-c",  NULL, ffOPTRD }
  };
#define NFILE asize(fnm)
  static int natoms=648;
  static bool bSymm=TRUE;
  static t_pargs pa[] = {
    { "-nat",  FALSE, etINT, { &natoms }, "Number of atoms" },
    { "-symm", FALSE, etBOOL,{ &bSymm  }, "Symmetrize the matrices" },
  };

  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
		    asize(desc),desc,0,NULL);

  bConf = (opt2bSet("-c",NFILE,fnm));
  if (bConf) {
    get_stx_coordnum (opt2fn("-c",NFILE,fnm),&natoms);
    init_t_atoms(&atoms,natoms,FALSE);
    snew(x,natoms);
    read_stx_conf(opt2fn("-c",NFILE,fnm),title,&atoms,x,NULL,box);
    set_pbc(&pbc,box);
  }
  snew(fn,2);
  fn[0] = opt2fn("-f1",NFILE,fnm);
  fn[1] = opt2fn("-f2",NFILE,fnm);
  
  snew(mat,2);  
  out = gmx_fio_fopen(ftp2fn(efOUT,NFILE,fnm),"w");
  mb  = sizeof(int)*sqr(natoms/1024.0);
  for(i=0; (i<2); i++) {
    in = gmx_fio_fopen(fn[i],"r");
    fprintf(stderr,"Reading %s\n",fn[i]);
    fprintf(out,   "Reading %s\n",fn[i]);
    fprintf(stderr,"Going to allocate %.0f Mb of memory\n",mb);
    fprintf(out,   "Going to allocate %.0f Mb of memory\n",mb);
    snew(mat[i],natoms);
    for(j=0; (j<natoms); j++) 
      snew(mat[i][j],natoms);
    nnb = read_nblist(in,out,mat[i],natoms,bSymm);
    gmx_fio_fclose(in);
    fprintf(stderr,"Interaction matrix %d has %d entries\n",i,nnb);
    fprintf(out,   "Interaction matrix %d has %d entries\n",i,nnb);
  }
  fprintf(stderr,"Comparing Interaction Matrices\n");
  mod=1;
  nmiss = 0;
  for(i=0; (i<natoms); i+=mod) {
    for(j=0; (j<natoms); j+=mod) {
      if (mat[0][i][j] != mat[1][i][j]) {
	fprintf(out,"i: %5d, j: %5d, shift[%s]: %3d, shift[%s]: %3d",
		i,j,fn[0],mat[0][i][j]-1,fn[1],mat[1][i][j]-1);
	if (bConf) {
	  pbc_dx(&pbc,x[i],x[j],dx);
	  fprintf(out," dist: %8.3f\n",norm(dx));
	}
	else
	  fprintf(out,"\n");
	nmiss++;
      }
    }
  }
  fprintf(out,"There were %d mismatches\n",nmiss);
  fprintf(out,"Done.\n");
  gmx_fio_fclose(out);
  fprintf(stderr,"There were %d mismatches\n",nmiss);
  fprintf(stderr,"Finished\n");
  
  thanx(stdout);
  
  return 0;
}

//.........这里部分代码省略.........
            a1  = molindex[mol+1];
            for (i = 0; i < nrefgrp; i++)
            {
                pbc_dx(&pbc, x[a0+srefat], bCom ? xref : x[index[0][i]], trial);
                rtry2 = norm2(trial);
                if (i == 0 || rtry2 < rdx2)
                {
                    copy_rvec(trial, dx);
                    rdx2 = rtry2;
                }
            }
            if (rdx2 < rcut2)
            {
                hist[static_cast<int>(std::sqrt(rdx2)*invbw)+1]++;
            }
            if (rdx2 >= rmin2 && rdx2 < rmax2)
            {
                unitv(dx, dx);
                clear_rvec(dip);
                qav = 0;
                for (a = a0; a < a1; a++)
                {
                    qav += atom[a].q;
                }
                qav /= (a1 - a0);
                for (a = a0; a < a1; a++)
                {
                    q = atom[a].q - qav;
                    for (d = 0; d < DIM; d++)
                    {
                        dip[d] += q*x[a][d];
                    }
                }
                for (d = 0; d < DIM; d++)
                {
                    dir[d] = -x[a0][d];
                }
                for (a = a0+1; a < a0+3; a++)
                {
                    for (d = 0; d < DIM; d++)
                    {
                        dir[d] += 0.5*x[a][d];
                    }
                }
                unitv(dir, dir);

                svmul(ENM2DEBYE, dip, dip);
                dip2   = norm2(dip);
                sdip  += std::sqrt(dip2);
                sdip2 += dip2;
                for (d = 0; d < DIM; d++)
                {
                    sinp  += dx[d]*dip[d];
                    sdinp += dx[d]*(dip[d] - refdip*dir[d]);
                }

                ntot++;
            }
        }
        nf++;

    }
    while (read_next_x(oenv, status, &t, x, box));

    gmx_rmpbc_done(gpbc);

    /* clean up */
    sfree(x);
    close_trj(status);

    fprintf(stderr, "Average number of molecules within %g nm is %.1f\n",
            rmax, static_cast<real>(ntot)/nf);
    if (ntot > 0)
    {
        sdip  /= ntot;
        sdip2 /= ntot;
        sinp  /= ntot;
        sdinp /= ntot;
        fprintf(stderr, "Average dipole:                               %f (D), std.dev. %f\n",
                sdip, std::sqrt(sdip2-sqr(sdip)));
        fprintf(stderr, "Average radial component of the dipole:       %f (D)\n",
                sinp);
        fprintf(stderr, "Average radial component of the polarization: %f (D)\n",
                sdinp);
    }

    fp = xvgropen(opt2fn("-o", NFILE, fnm),
                  "Cumulative solvent distribution", "r (nm)", "molecules", oenv);
    nmol = 0;
    for (i = 0; i <= nbin; i++)
    {
        nmol += hist[i];
        fprintf(fp, "%g %g\n", i*bw, nmol/nf);
    }
    xvgrclose(fp);

    do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);

    return 0;
}

//.........这里部分代码省略.........
    { "-seed",   FALSE, etINT , {&seed},      
      "random generator seed"},
    { "-vdwd",   FALSE, etREAL, {&r_distance},
      "default vdwaals distance"},
    { "-shell",  FALSE, etREAL, {&r_shell},
      "thickness of optional water layer around solute" },
    { "-maxsol", FALSE, etINT,  {&max_sol},
      "maximum number of solvent molecules to add if they fit in the box. If zero (default) this is ignored" },
    { "-vel",    FALSE, etBOOL, {&bReadV},
      "keep velocities from input solute and solvent" }
  };

  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv, PCA_BE_NICE,NFILE,fnm,asize(pa),pa,
		    asize(desc),desc,asize(bugs),bugs,&oenv);
  
  bInsert   = opt2bSet("-ci",NFILE,fnm) && (nmol_ins > 0);
  bSol      = opt2bSet("-cs",NFILE,fnm);
  bProt     = opt2bSet("-cp",NFILE,fnm);
  bBox      = opt2parg_bSet("-box",asize(pa),pa);
     
  /* check input */
  if (bInsert && nmol_ins<=0)
    gmx_fatal(FARGS,"When specifying inserted molecules (-ci), "
		"-nmol must be larger than 0");
  if (!bProt && !bBox)
    gmx_fatal(FARGS,"When no solute (-cp) is specified, "
		"a box size (-box) must be specified");

  aps = gmx_atomprop_init();
  
  if (bProt) {
    /*generate a solute configuration */
    conf_prot = opt2fn("-cp",NFILE,fnm);
    title = read_prot(conf_prot,&atoms,&x,bReadV?&v:NULL,&r,&ePBC,box,
		      aps,r_distance);
    if (bReadV && !v)
      fprintf(stderr,"Note: no velocities found\n");
    if (atoms.nr == 0) {
      fprintf(stderr,"Note: no atoms in %s\n",conf_prot);
      bProt = FALSE;
    }
  } 
  if (!bProt) {
    atoms.nr=0;
    atoms.nres=0;
    atoms.resinfo=NULL;
    atoms.atomname=NULL;
    atoms.atom=NULL;
    atoms.pdbinfo=NULL;
    x=NULL;
    r=NULL;
  }
  if (bBox) {
    ePBC = epbcXYZ;
    clear_mat(box);
    box[XX][XX]=new_box[XX];
    box[YY][YY]=new_box[YY];
    box[ZZ][ZZ]=new_box[ZZ];
  }
  if (det(box) == 0) 
    gmx_fatal(FARGS,"Undefined solute box.\nCreate one with editconf "
		"or give explicit -box command line option");
  
  /* add nmol_ins molecules of atoms_ins 
     in random orientation at random place */

//.........这里部分代码省略.........
    { "-oneacf", FALSE, etBOOL, {&bAverCorr},
      "Calculate one ACF over all sets" }
  };
#define NPA asize(pa)

  FILE     *out,*out_fit;
  int      n,nlast,s,nset,i,j=0;
  real     **val,*t,dt,tot,error;
  double   *av,*sig,cum1,cum2,cum3,cum4,db;
  char     *acfile,*msdfile,*ccfile,*distfile,*avfile,*eefile,*fitfile;
  
  t_filenm fnm[] = { 
    { efXVG, "-f",    "graph",    ffREAD   },
    { efXVG, "-ac",   "autocorr", ffOPTWR  },
    { efXVG, "-msd",  "msd",      ffOPTWR  },
    { efXVG, "-cc",   "coscont",  ffOPTWR  },
    { efXVG, "-dist", "distr",    ffOPTWR  },
    { efXVG, "-av",   "average",  ffOPTWR  },
    { efXVG, "-ee",   "errest",   ffOPTWR  },
    { efLOG, "-g",    "fitlog",   ffOPTWR  }
  }; 
#define NFILE asize(fnm) 

  int     npargs;
  t_pargs *ppa;

  npargs = asize(pa); 
  ppa    = add_acf_pargs(&npargs,pa);
  
  CopyRight(stderr,argv[0]); 
  parse_common_args(&argc,argv,PCA_CAN_VIEW,
		    NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL); 

  acfile   = opt2fn_null("-ac",NFILE,fnm);
  msdfile  = opt2fn_null("-msd",NFILE,fnm);
  ccfile   = opt2fn_null("-cc",NFILE,fnm);
  distfile = opt2fn_null("-dist",NFILE,fnm);
  avfile   = opt2fn_null("-av",NFILE,fnm);
  eefile   = opt2fn_null("-ee",NFILE,fnm);
  if (opt2parg_bSet("-fitfn",npargs,ppa)) 
    fitfile  = opt2fn("-g",NFILE,fnm);
  else
    fitfile  = opt2fn_null("-g",NFILE,fnm);
    
  val=read_xvg_time(opt2fn("-f",NFILE,fnm),bHaveT,
		    opt2parg_bSet("-b",npargs,ppa),tb,
		    opt2parg_bSet("-e",npargs,ppa),te,
		    nsets_in,&nset,&n,&dt,&t);
  printf("Read %d sets of %d points, dt = %g\n\n",nset,n,dt);
  
  if (bDer) {
    printf("Calculating the derivative as (f[i+%d]-f[i])/(%d*dt)\n\n",
	    d,d);
    n -= d;
    for(s=0; s<nset; s++)
      for(i=0; (i<n); i++)
	val[s][i] = (val[s][i+d]-val[s][i])/(d*dt);
  }
  if (bIntegrate) {
    real sum,stddev;
    printf("Calculating the integral using the trapezium rule\n");
    
    if (bXYdy) {
      sum = evaluate_integral(n,t,val[0],val[1],aver_start,&stddev);
      printf("Integral %10.3f +/- %10.5f\n",sum,stddev);
    }

//.........这里部分代码省略.........
    int             n, e_n, nlast, s, nset, e_nset, d_nset, i, j = 0, *idim, *ibox;
    real          **val, **et_val, *t, *e_t, e_dt, d_dt, dt, tot, error;
    real           *rmin, *rmax;
    double         *av, *sig, cum1, cum2, cum3, cum4, db;
    const char     *fn_ge, *fn_ene;
    output_env_t    oenv;
    gmx_int64_t     num_grid_points;

    t_filenm        fnm[] = {
        { efXVG, "-f",    "graph",    ffREAD   },
        { efXVG, "-ge",   "gibbs",    ffOPTRD  },
        { efXVG, "-ene",  "esham",    ffOPTRD  },
        { efXVG, "-dist", "ener",     ffOPTWR  },
        { efXVG, "-histo", "edist",    ffOPTWR  },
        { efNDX, "-bin",  "bindex",   ffOPTWR  },
        { efXPM, "-lp",   "prob",     ffOPTWR  },
        { efXPM, "-ls",   "gibbs",    ffOPTWR  },
        { efXPM, "-lsh",  "enthalpy", ffOPTWR  },
        { efXPM, "-lss",  "entropy",  ffOPTWR  },
        { efPDB, "-ls3",  "gibbs3",   ffOPTWR  },
        { efLOG, "-g",    "shamlog",  ffOPTWR  }
    };
#define NFILE asize(fnm)

    int     npargs;

    npargs = asize(pa);
    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
                           NFILE, fnm, npargs, pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }

    val = read_xvg_time(opt2fn("-f", NFILE, fnm), bHaveT,
                        opt2parg_bSet("-b", npargs, pa), tb-ttol,
                        opt2parg_bSet("-e", npargs, pa), te+ttol,
                        nsets_in, &nset, &n, &dt, &t);
    printf("Read %d sets of %d points, dt = %g\n\n", nset, n, dt);

    fn_ge  = opt2fn_null("-ge", NFILE, fnm);
    fn_ene = opt2fn_null("-ene", NFILE, fnm);

    if (fn_ge && fn_ene)
    {
        gmx_fatal(FARGS, "Can not do free energy and energy corrections at the same time");
    }

    if (fn_ge || fn_ene)
    {
        et_val = read_xvg_time(fn_ge ? fn_ge : fn_ene, bHaveT,
                               opt2parg_bSet("-b", npargs, pa), tb-ttol,
                               opt2parg_bSet("-e", npargs, pa), te+ttol,
                               1, &e_nset, &e_n, &e_dt, &e_t);
        if (fn_ge)
        {
            if (e_nset != 1)
            {
                gmx_fatal(FARGS, "Can only handle one free energy component in %s",
                          fn_ge);
            }
        }
        else
        {
            if (e_nset != 1 && e_nset != 2)
            {
                gmx_fatal(FARGS, "Can only handle one energy component or one energy and one T in %s",

int gmx_nmens(int argc,char *argv[])
{
  static char *desc[] = {
    "[TT]g_nmens[tt] generates an ensemble around an average structure",
    "in a subspace which is defined by a set of normal modes (eigenvectors).",
    "The eigenvectors are assumed to be mass-weighted.",
    "The position along each eigenvector is randomly taken from a Gaussian",
    "distribution with variance kT/eigenvalue.[PAR]",
    "By default the starting eigenvector is set to 7, since the first six",
    "normal modes are the translational and rotational degrees of freedom." 
  };
  static int  nstruct=100,first=7,last=-1,seed=-1;
  static real temp=300.0;
  t_pargs pa[] = {
    { "-temp",  FALSE, etREAL, {&temp}, 
      "Temperature in Kelvin" },
    { "-seed", FALSE, etINT, {&seed},     
      "Random seed, -1 generates a seed from time and pid" },
    { "-num", FALSE, etINT, {&nstruct},     
      "Number of structures to generate" },
    { "-first", FALSE, etINT, {&first},     
      "First eigenvector to use (-1 is select)" },
    { "-last",  FALSE, etINT, {&last}, 
      "Last eigenvector to use (-1 is till the last)" }
  };
#define NPA asize(pa)
  
  int        out;
  int        status,trjout;
  t_topology top;
  int        ePBC;
  t_atoms    *atoms;
  rvec       *xtop,*xref,*xav,*xout1,*xout2;
  bool       bDMR,bDMA,bFit;
  int        nvec,*eignr=NULL;
  rvec       **eigvec=NULL;
  matrix     box;
  real       *eigval,totmass,*invsqrtm,t,disp;
  int        natoms,neigval;
  char       *grpname,*indexfile,title[STRLEN];
  int        i,j,d,s,v;
  int        nout,*iout,noutvec,*outvec;
  atom_id    *index;
  real       rfac,invfr,rhalf,jr;
  int *      eigvalnr;
  
  unsigned long      jran;
  const unsigned long im = 0xffff;
  const unsigned long ia = 1093;
  const unsigned long ic = 18257;

  t_filenm fnm[] = { 
    { efTRN, "-v",    "eigenvec",    ffREAD  },
    { efXVG, "-e",    "eigenval",    ffREAD  },
    { efTPS, NULL,    NULL,          ffREAD },
    { efNDX, NULL,    NULL,          ffOPTRD },
    { efTRO, "-o",    "ensemble",    ffWRITE }
  }; 
#define NFILE asize(fnm) 

  CopyRight(stderr,argv[0]); 
  parse_common_args(&argc,argv,PCA_BE_NICE,
		    NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL); 

  indexfile=ftp2fn_null(efNDX,NFILE,fnm);

  read_eigenvectors(opt2fn("-v",NFILE,fnm),&natoms,&bFit,
		    &xref,&bDMR,&xav,&bDMA,&nvec,&eignr,&eigvec,&eigval);

  read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xtop,NULL,box,bDMA);
  atoms=&top.atoms;

  printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n",natoms);
  get_index(atoms,indexfile,1,&i,&index,&grpname);
  if (i!=natoms)
    gmx_fatal(FARGS,"you selected a group with %d elements instead of %d",
		i,natoms);
  printf("\n");
  
  snew(invsqrtm,natoms);
  if (bDMA) {
    for(i=0; (i<natoms); i++)
      invsqrtm[i] = invsqrt(atoms->atom[index[i]].m);
  } else {
    for(i=0; (i<natoms); i++)
      invsqrtm[i]=1.0;
  }
  
  if (last==-1)
      last=natoms*DIM;
  if (first>-1) 
  {
      /* make an index from first to last */
      nout=last-first+1;
      snew(iout,nout);
      for(i=0; i<nout; i++)
          iout[i]=first-1+i;
  }
  else 
  {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                " calculating fractions as defined in this program\n"
                "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
        bFrac = FALSE;
    }

    if ( (bTrans || bFrac || bCorr) && mult == 3)
    {
        gmx_fatal(FARGS, "Can only do transition, fraction or correlation\n"
                  "on dihedrals. Select -d\n");
    }

    /*
     * We need to know the nr of frames so we can allocate memory for an array
     * with all dihedral angles at all timesteps. Works for me.
     */
    if (bTrans || bCorr  || bALL || opt2bSet("-or", NFILE, fnm))
    {
        snew(dih, nangles);
    }

    snew(angstat, maxangstat);

    read_ang_dih(ftp2fn(efTRX, NFILE, fnm), (mult == 3),
                 bALL || bCorr || bTrans || opt2bSet("-or", NFILE, fnm),
                 bRb, bPBC, maxangstat, angstat,
                 &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
                 oenv);

    dt = (time[nframes-1]-time[0])/(nframes-1);

    if (bAver)
    {
        sprintf(title, "Average Angle: %s", grpname);
        out = xvgropen(opt2fn("-ov", NFILE, fnm),
                       title, "Time (ps)", "Angle (degrees)", oenv);
        for (i = 0; (i < nframes); i++)
        {
            fprintf(out, "%10.5f  %8.3f", time[i], aver_angle[i]*RAD2DEG);
            if (bALL)
            {
                for (j = 0; (j < nangles); j++)
                {
                    if (bPBC)
                    {
                        real dd = dih[j][i];
                        fprintf(out, "  %8.3f", std::atan2(std::sin(dd), std::cos(dd))*RAD2DEG);
                    }
                    else
                    {
                        fprintf(out, "  %8.3f", dih[j][i]*RAD2DEG);
                    }
                }
            }
            fprintf(out, "\n");
        }
        xvgrclose(out);
    }
    if (opt2bSet("-or", NFILE, fnm))
    {
        dump_dih_trr(nframes, nangles, dih, opt2fn("-or", NFILE, fnm), time);
    }

    if (bFrac)
    {
        sprintf(title, "Trans fraction: %s", grpname);
        out = xvgropen(opt2fn("-of", NFILE, fnm),

int gmx_dyndom(int argc, char *argv[])
{
    const char  *desc[] = {
        "[THISMODULE] reads a [REF].pdb[ref] file output from DynDom",
        "(http://www.cmp.uea.ac.uk/dyndom/).",
        "It reads the coordinates, the coordinates of the rotation axis,",
        "and an index file containing the domains.",
        "Furthermore, it takes the first and last atom of the arrow file",
        "as command line arguments (head and tail) and",
        "finally it takes the translation vector (given in DynDom info file)",
        "and the angle of rotation (also as command line arguments). If the angle",
        "determined by DynDom is given, one should be able to recover the",
        "second structure used for generating the DynDom output.",
        "Because of limited numerical accuracy this should be verified by",
        "computing an all-atom RMSD (using [gmx-confrms]) rather than by file",
        "comparison (using diff).[PAR]",
        "The purpose of this program is to interpolate and extrapolate the",
        "rotation as found by DynDom. As a result unphysical structures with",
        "long or short bonds, or overlapping atoms may be produced. Visual",
        "inspection, and energy minimization may be necessary to",
        "validate the structure."
    };
    static real  trans0 = 0;
    static rvec  head   = { 0, 0, 0 };
    static rvec  tail   = { 0, 0, 0 };
    static real  angle0 = 0, angle1 = 0, maxangle = 0;
    static int   label  = 0, nframes = 11;
    t_pargs      pa[]   = {
        { "-firstangle",    FALSE, etREAL, {&angle0},
          "Angle of rotation about rotation vector" },
        { "-lastangle",    FALSE, etREAL, {&angle1},
          "Angle of rotation about rotation vector" },
        { "-nframe",   FALSE, etINT,  {&nframes},
          "Number of steps on the pathway" },
        { "-maxangle", FALSE, etREAL, {&maxangle},
          "DymDom dtermined angle of rotation about rotation vector" },
        { "-trans",    FALSE, etREAL, {&trans0},
          "Translation (Angstrom) along rotation vector (see DynDom info file)" },
        { "-head",     FALSE, etRVEC, {head},
          "First atom of the arrow vector" },
        { "-tail",     FALSE, etRVEC, {tail},
          "Last atom of the arrow vector" }
    };
    int          i, j, natoms, isize;
    t_trxstatus *status;
    atom_id     *index = NULL, *index_all;
    char         title[256], *grpname;
    t_atoms      atoms;
    real         angle, trans;
    rvec        *x, *v, *xout, *vout;
    matrix       box;
    output_env_t oenv;

    t_filenm     fnm[] = {
        { efPDB, "-f", "dyndom",  ffREAD },
        { efTRO, "-o", "rotated", ffWRITE },
        { efNDX, "-n", "domains", ffREAD }
    };
#define NFILE asize(fnm)

    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
                           asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }

    if (maxangle == 0)
    {
        gmx_fatal(FARGS, "maxangle not given");
    }

    get_stx_coordnum (opt2fn("-f", NFILE, fnm), &natoms);
    init_t_atoms(&atoms, natoms, TRUE);
    snew(x, natoms);
    snew(v, natoms);
    read_stx_conf(opt2fn("-f", NFILE, fnm), title, &atoms, x, v, NULL, box);
    snew(xout, natoms);
    snew(vout, natoms);

    printf("Select group to rotate:\n");
    rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
    printf("Going to rotate %s containg %d atoms\n", grpname, isize);

    snew(index_all, atoms.nr);
    for (i = 0; (i < atoms.nr); i++)
    {
        index_all[i] = i;
    }

    status = open_trx(opt2fn("-o", NFILE, fnm), "w");

    label = 'A';
    for (i = 0; (i < nframes); i++, label++)
    {
        angle = angle0 + (i*(angle1-angle0))/(nframes-1);
        trans = trans0*0.1*angle/maxangle;
        printf("Frame: %2d (label %c), angle: %8.3f deg., trans: %8.3f nm\n",
               i, label, angle, trans);
        rot_conf(&atoms, x, v, trans, angle, head, tail, isize, index, xout, vout);

//.........这里部分代码省略.........

int gmx_morph(int argc, char *argv[])
{
    const char      *desc[] = {
        "[THISMODULE] does a linear interpolation of conformations in order to",
        "create intermediates. Of course these are completely unphysical, but",
        "that you may try to justify yourself. Output is in the form of a ",
        "generic trajectory. The number of intermediates can be controlled with",
        "the [TT]-ninterm[tt] flag. The first and last flag correspond to the way of",
        "interpolating: 0 corresponds to input structure 1 while",
        "1 corresponds to input structure 2.",
        "If you specify [TT]-first[tt] < 0 or [TT]-last[tt] > 1 extrapolation will be",
        "on the path from input structure x[SUB]1[sub] to x[SUB]2[sub]. In general, the coordinates",
        "of the intermediate x(i) out of N total intermediates correspond to:[PAR]",
        "x(i) = x[SUB]1[sub] + (first+(i/(N-1))*(last-first))*(x[SUB]2[sub]-x[SUB]1[sub])[PAR]",
        "Finally the RMSD with respect to both input structures can be computed",
        "if explicitly selected ([TT]-or[tt] option). In that case, an index file may be",
        "read to select the group from which the RMS is computed."
    };
    t_filenm         fnm[] = {
        { efSTX, "-f1", "conf1",  ffREAD },
        { efSTX, "-f2", "conf2",  ffREAD },
        { efTRX, "-o",  "interm", ffWRITE },
        { efXVG, "-or", "rms-interm", ffOPTWR },
        { efNDX, "-n",  "index",  ffOPTRD }
    };
#define NFILE asize(fnm)
    static  int      ninterm = 11;
    static  real     first   = 0.0;
    static  real     last    = 1.0;
    static  gmx_bool bFit    = TRUE;
    t_pargs          pa []   = {
        { "-ninterm", FALSE, etINT,  {&ninterm},
          "Number of intermediates" },
        { "-first",   FALSE, etREAL, {&first},
          "Corresponds to first generated structure (0 is input x[SUB]1[sub], see above)" },
        { "-last",    FALSE, etREAL, {&last},
          "Corresponds to last generated structure (1 is input x[SUB]2[sub], see above)" },
        { "-fit",     FALSE, etBOOL, {&bFit},
          "Do a least squares fit of the second to the first structure before interpolating" }
    };
    const char      *leg[] = { "Ref = 1\\Sst\\N conf", "Ref = 2\\Snd\\N conf" };
    FILE            *fp    = NULL;
    int              i, isize, is_lsq, nat1, nat2;
    t_trxstatus     *status;
    atom_id         *index, *index_lsq, *index_all, *dummy;
    t_atoms          atoms;
    rvec            *x1, *x2, *xx, *v;
    matrix           box;
    real             rms1, rms2, fac, *mass;
    char             title[STRLEN], *grpname;
    gmx_bool         bRMS;
    output_env_t     oenv;

    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc,
                           0, NULL, &oenv))
    {
        return 0;
    }
    get_stx_coordnum (opt2fn("-f1", NFILE, fnm), &nat1);
    get_stx_coordnum (opt2fn("-f2", NFILE, fnm), &nat2);
    if (nat1 != nat2)
    {
        gmx_fatal(FARGS, "Number of atoms in first structure is %d, in second %d",
                  nat1, nat2);
    }

    init_t_atoms(&atoms, nat1, TRUE);
    snew(x1, nat1);
    snew(x2, nat1);
    snew(xx, nat1);
    snew(v, nat1);

    read_stx_conf(opt2fn("-f1", NFILE, fnm), title, &atoms, x1, v, NULL, box);
    read_stx_conf(opt2fn("-f2", NFILE, fnm), title, &atoms, x2, v, NULL, box);

    snew(mass, nat1);
    snew(index_all, nat1);
    for (i = 0; (i < nat1); i++)
    {
        mass[i]      = 1;
        index_all[i] = i;
    }
    if (bFit)
    {
        printf("Select group for LSQ superposition:\n");
        get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &is_lsq, &index_lsq,
                  &grpname);
        reset_x(is_lsq, index_lsq, nat1, index_all, x1, mass);
        reset_x(is_lsq, index_lsq, nat1, index_all, x2, mass);
        do_fit(nat1, mass, x1, x2);
    }

    bRMS = opt2bSet("-or", NFILE, fnm);
    if (bRMS)
    {
        fp = xvgropen(opt2fn("-or", NFILE, fnm), "RMSD", "Conf", "(nm)", oenv);
        xvgr_legend(fp, asize(leg), leg, oenv);
        printf("Select group for RMSD calculation:\n");
        get_index(&atoms, opt2fn_null("-n", NFILE, fnm), 1, &isize, &index, &grpname);
//.........这里部分代码省略.........

int gmx_clustsize(int argc,char *argv[])
{
  static char *desc[] = {
    "This program computes the size distributions of molecular/atomic clusters in",
    "the gas phase. The output is given in the form of a XPM file.",
    "The total number of clusters is written to a XVG file.[PAR]",
    "When the [TT]-mol[tt] option is given clusters will be made out of",
    "molecules rather than atoms, which allows clustering of large molecules.",
    "In this case an index file would still contain atom numbers",
    "or your calculcation will die with a SEGV.[PAR]",
    "When velocities are present in your trajectory, the temperature of",
    "the largest cluster will be printed in a separate xvg file assuming",
    "that the particles are free to move. If you are using constraints,",
    "please correct the temperature. For instance water simulated with SHAKE",
    "or SETTLE will yield a temperature that is 1.5 times too low. You can",
    "compensate for this with the -ndf option. Remember to take the removal",
    "of center of mass motion into account.[PAR]",
    "The [TT]-mc[tt] option will produce an index file containing the",
    "atom numbers of the largest cluster."
  };
  
  static real cutoff   = 0.35;
  static int  nskip    = 0;
  static int  nlevels  = 20;
  static int  ndf      = -1;
  static bool bMol     = FALSE;
  static bool bPBC     = TRUE;
  static rvec rlo      = { 1.0, 1.0, 0.0 };
  static rvec rhi      = { 0.0, 0.0, 1.0 };
  t_pargs pa[] = {
    { "-cut",      FALSE, etREAL, {&cutoff},
      "Largest distance (nm) to be considered in a cluster" },
    { "-mol",      FALSE, etBOOL, {&bMol},
      "Cluster molecules rather than atoms (needs tpr file)" },
    { "-pbc",      FALSE, etBOOL, {&bPBC},
      "Use periodic boundary conditions" },
    { "-nskip",    FALSE, etINT,  {&nskip},
      "Number of frames to skip between writing" },
    { "-nlevels",  FALSE, etINT,  {&nlevels},
      "Number of levels of grey in xpm output" },
    { "-ndf",      FALSE, etINT,  {&ndf},
      "Number of degrees of freedom of the entire system for temperature calculation. If not set the number of atoms times three is used." },
    { "-rgblo",    FALSE, etRVEC, {rlo},
      "RGB values for the color of the lowest occupied cluster size" },
    { "-rgbhi",    FALSE, etRVEC, {rhi},
      "RGB values for the color of the highest occupied cluster size" }
  };
#define NPA asize(pa)
  char       *fnNDX,*fnTPR;
  bool       bSQ,bRDF;
  t_rgb      rgblo,rgbhi;
  
  t_filenm   fnm[] = {
    { efTRX, "-f",  NULL,         ffREAD  },
    { efTPR, NULL,  NULL,         ffOPTRD },
    { efNDX, NULL,  NULL,         ffOPTRD },
    { efXPM, "-o", "csize",       ffWRITE },
    { efXPM, "-ow","csizew",      ffWRITE },
    { efXVG, "-nc","nclust",      ffWRITE },
    { efXVG, "-mc","maxclust",    ffWRITE },
    { efXVG, "-ac","avclust",     ffWRITE },
    { efXVG, "-hc","histo-clust", ffWRITE },
    { efXVG, "-temp","temp",     ffOPTWR },
    { efNDX, "-mcn", "maxclust", ffOPTWR }
  };
#define NFILE asize(fnm)
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
		    NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);

  fnNDX = ftp2fn_null(efNDX,NFILE,fnm);
  rgblo.r = rlo[XX],rgblo.g = rlo[YY],rgblo.b = rlo[ZZ];
  rgbhi.r = rhi[XX],rgbhi.g = rhi[YY],rgbhi.b = rhi[ZZ];
    
  fnTPR = ftp2fn_null(efTPR,NFILE,fnm);
  if (bMol && !fnTPR)
    gmx_fatal(FARGS,"You need a tpr file for the -mol option");
    
  clust_size(fnNDX,ftp2fn(efTRX,NFILE,fnm),opt2fn("-o",NFILE,fnm),
	     opt2fn("-ow",NFILE,fnm),
	     opt2fn("-nc",NFILE,fnm),opt2fn("-ac",NFILE,fnm),
	     opt2fn("-mc",NFILE,fnm),opt2fn("-hc",NFILE,fnm),
	     opt2fn("-temp",NFILE,fnm),opt2fn("-mcn",NFILE,fnm),
	     bMol,bPBC,fnTPR,
	     cutoff,nskip,nlevels,rgblo,rgbhi,ndf);

  thanx(stderr);
  
  return 0;
}

int gmx_principal(int argc, char *argv[])
{
    const char     *desc[] = {
        "[TT]g_principal[tt] calculates the three principal axes of inertia for a group",
        "of atoms.",
    };
    static gmx_bool foo = FALSE;

    t_pargs         pa[] = {
        { "-foo",      FALSE, etBOOL, {&foo}, "Dummy option to avoid empty array" }
    };
    t_trxstatus    *status;
    t_topology      top;
    int             ePBC;
    real            t;
    rvec      *     x;

    int             natoms;
    char           *grpname, title[256];
    int             i, j, m, gnx, nam, mol;
    atom_id        *index;
    rvec            a1, a2, a3, moi;
    FILE      *     axis1;
    FILE      *     axis2;
    FILE      *     axis3;
    FILE      *     fmoi;
    matrix          axes, box;
    output_env_t    oenv;
    gmx_rmpbc_t     gpbc = NULL;


    t_filenm fnm[] = {
        { efTRX, "-f",   NULL,       ffREAD },
        { efTPS, NULL,   NULL,       ffREAD },
        { efNDX, NULL,   NULL,       ffOPTRD },
        { efDAT, "-a1",  "axis1",    ffWRITE },
        { efDAT, "-a2",  "axis2",    ffWRITE },
        { efDAT, "-a3",  "axis3",    ffWRITE },
        { efDAT, "-om",  "moi",      ffWRITE }
    };
#define NFILE asize(fnm)

    parse_common_args(&argc, argv,
                      PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
                      NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);

    axis1 = ffopen(opt2fn("-a1", NFILE, fnm), "w");
    axis2 = ffopen(opt2fn("-a2", NFILE, fnm), "w");
    axis3 = ffopen(opt2fn("-a3", NFILE, fnm), "w");
    fmoi  = ffopen(opt2fn("-om", NFILE, fnm), "w");

    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box, TRUE);

    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);

    natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);

    gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);

    do
    {
        gmx_rmpbc(gpbc, natoms, box, x);

        calc_principal_axes(&top, x, index, gnx, axes, moi);

        fprintf(axis1, "%15.10f     %15.10f  %15.10f  %15.10f\n", t, axes[XX][XX], axes[YY][XX], axes[ZZ][XX]);
        fprintf(axis2, "%15.10f     %15.10f  %15.10f  %15.10f\n", t, axes[XX][YY], axes[YY][YY], axes[ZZ][YY]);
        fprintf(axis3, "%15.10f     %15.10f  %15.10f  %15.10f\n", t, axes[XX][ZZ], axes[YY][ZZ], axes[ZZ][ZZ]);
        fprintf(fmoi,  "%15.10f     %15.10f  %15.10f  %15.10f\n", t, moi[XX], moi[YY], moi[ZZ]);
    }
    while (read_next_x(oenv, status, &t, x, box));

    gmx_rmpbc_done(gpbc);


    close_trj(status);
    ffclose(axis1);
    ffclose(axis2);
    ffclose(axis3);
    ffclose(fmoi);

    return 0;
}