C++ please_cite函数代码示例

int main(int argc,char *argv[])
{
  const char *desc[] = {
    "[TT]g_anadock[tt] analyses the results of an Autodock run and clusters the",
    "structures together, based on distance or RMSD. The docked energy",
    "and free energy estimates are analysed, and for each cluster the",
    "energy statistics are printed.[PAR]",
    "An alternative approach to this is to cluster the structures first",
    "using [TT]g_cluster[tt] and then sort the clusters on either lowest",
    "energy or average energy."
  };
  t_filenm fnm[] = {
    { efPDB, "-f", NULL,       ffREAD  },
    { efPDB, "-ox", "cluster", ffWRITE },
    { efXVG, "-od", "edocked", ffWRITE },
    { efXVG, "-of", "efree",   ffWRITE },
    { efLOG, "-g",  "anadock", ffWRITE }
  };
  output_env_t oenv;
#define NFILE asize(fnm)
  static gmx_bool bFree=FALSE,bRMS=TRUE;
  static real cutoff = 0.2;
  t_pargs pa[] = {
    { "-free",   FALSE, etBOOL, {&bFree}, 
      "Use Free energy estimate from autodock for sorting the classes" },
    { "-rms",    FALSE, etBOOL, {&bRMS},
      "Cluster on RMS or distance" },
    { "-cutoff", FALSE, etREAL, {&cutoff},
      "Maximum RMSD/distance for belonging to the same cluster" }
  };
#define NPA asize(pa)

  FILE      *fp;
  t_pdbfile **pdbf=NULL;
  int       npdbf;
  
  CopyRight(stderr,argv[0]);
  parse_common_args(&argc,argv,0,NFILE,fnm,NPA,pa, asize(desc),desc,0,
                    NULL,&oenv);
  
  fp = ffopen(opt2fn("-g",NFILE,fnm),"w");
  please_cite(stdout,"Hetenyi2002b");
  please_cite(fp,"Hetenyi2002b");
  
  pdbf = read_em_all(opt2fn("-f",NFILE,fnm),&npdbf);
  
  analyse_em_all(npdbf,pdbf,opt2fn("-od",NFILE,fnm),opt2fn("-of",NFILE,fnm),
                 oenv);
  
  cluster_em_all(fp,npdbf,pdbf,opt2fn("-ox",NFILE,fnm),
		 bFree,bRMS,cutoff);
  
  thanx(fp);
  ffclose(fp);
  
  thanx(stdout);
  
  return 0;
}

static void vdw_warning(FILE *fp)
{
    if (NULL != fp)
    {
        fprintf(fp, "NOTE: From version 5.0 %s uses the Van der Waals radii\n",
                gmx::getProgramContext().displayName());
        fprintf(fp, "from the source below. This means the results may be different\n");
        fprintf(fp, "compared to previous GROMACS versions.\n");
        please_cite(fp, "Bondi1964a");
    }
}

void calc_rffac(FILE *fplog,int eel,real eps_r,real eps_rf,real Rc,real Temp,
		real zsq,matrix box,
		real *kappa,real *krf,real *crf)
{
  /* Compute constants for Generalized reaction field */
  real   k1,k2,I,vol,rmin;
  
  if (EEL_RF(eel)) {
    vol     = det(box);
    if (eel == eelGRF) {
      /* Consistency check */
      if (Temp <= 0.0)
	gmx_fatal(FARGS,"Temperature is %f while using"
		    " Generalized Reaction Field\n",Temp);
      /* Ionic strength (only needed for eelGRF */
      I       = 0.5*zsq/vol;
      *kappa  = sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
    }
    else {
      I      = 0;
      *kappa = 0;
    }

    /* eps == 0 signals infinite dielectric */
    if (eps_rf == 0) {
      *krf = 1/(2*Rc*Rc*Rc);
    } else {
      k1   = 1 + *kappa*Rc;
      k2   = eps_rf*sqr((real)(*kappa*Rc));
      
      *krf = ((eps_rf - eps_r)*k1 + 0.5*k2)/((2*eps_rf + eps_r)*k1 + k2)/(Rc*Rc*Rc);
    }
    *crf   = 1/Rc + *krf*Rc*Rc;
    rmin   = pow(*krf*2.0,-1.0/3.0);
    
    if (fplog) {
      if (eel == eelGRF) {
	please_cite(fplog,"Tironi95a");
	fprintf(fplog,"%s:\n"
		"epsRF = %10g, I   = %10g, volume = %10g, kappa  = %10g\n"
		"rc    = %10g, krf = %10g, crf    = %10g, epsfac = %10g\n",
		eel_names[eel],eps_rf,I,vol,*kappa,Rc,*krf,*crf,
		ONE_4PI_EPS0/eps_r);
      } else {
	fprintf(fplog,"%s:\n"
		"epsRF = %g, rc = %g, krf = %g, crf = %g, epsfac = %g\n",
		eel_names[eel],eps_rf,Rc,*krf,*crf,ONE_4PI_EPS0/eps_r);
      }
      fprintf(fplog,
	      "The electrostatics potential has its minimum at r = %g\n",
	      rmin);
    }
  }
}

int gmx_saxs(int argc, char *argv[])
{
    const char  *desc[] = {
        "[THISMODULE] calculates SAXS structure factors for given index",
        "groups based on Cromer's method.",
        "Both topology and trajectory files are required."
    };

    static real  start_q = 0.0, end_q = 60.0, energy = 12.0;
    static int   ngroups = 1;

    t_pargs      pa[] = {
        { "-ng",       FALSE, etINT, {&ngroups},
          "Number of groups to compute SAXS" },
        {"-startq", FALSE, etREAL, {&start_q},
         "Starting q (1/nm) "},
        {"-endq", FALSE, etREAL, {&end_q},
         "Ending q (1/nm)"},
        {"-energy", FALSE, etREAL, {&energy},
         "Energy of the incoming X-ray (keV) "}
    };
#define NPA asize(pa)
    const char  *fnTPS, *fnTRX, *fnNDX, *fnDAT = NULL;
    output_env_t oenv;

    t_filenm     fnm[] = {
        { efTRX, "-f",  NULL,      ffREAD },
        { efTPS, NULL,  NULL,      ffREAD },
        { efNDX, NULL,  NULL,      ffOPTRD },
        { efDAT, "-d",  "sfactor", ffOPTRD },
        { efXVG, "-sq", "sq",      ffWRITE },
    };
#define NFILE asize(fnm)
    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
                           NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }

    fnTPS = ftp2fn(efTPS, NFILE, fnm);
    fnTRX = ftp2fn(efTRX, NFILE, fnm);
    fnDAT = ftp2fn(efDAT, NFILE, fnm);
    fnNDX = ftp2fn_null(efNDX, NFILE, fnm);

    do_scattering_intensity(fnTPS, fnNDX, opt2fn("-sq", NFILE, fnm),
                            fnTRX, fnDAT,
                            start_q, end_q, energy, ngroups, oenv);

    please_cite(stdout, "Cromer1968a");

    return 0;
}

void do_pp2shifts(FILE *fp, int nf, int nlist, t_dlist dlist[], real **dih)
{
    t_shiftdata *ca_sd, *co_sd, *ha_sd, *cb_sd;
    int          i, j, Phi, Psi;
    real         phi, psi;
    real         ca, co, ha, cb;

    /* Read the shift files */
    ca_sd = read_shifts("ca-shift.dat");
    cb_sd = read_shifts("cb-shift.dat");
    ha_sd = read_shifts("ha-shift.dat");
    co_sd = read_shifts("co-shift.dat");

    fprintf(fp, "\n *** Chemical shifts from the chemical shift index ***\n");
    please_cite(fp, "Wishart98a");
    fprintf(fp, "%12s  %10s  %10s  %10s  %10s\n",
            "Residue", "delta Ca", "delta Ha", "delta CO", "delta Cb");
    for (i = 0; (i < nlist); i++)
    {
        if ((has_dihedral(edPhi, &(dlist[i]))) &&
            (has_dihedral(edPsi, &(dlist[i]))))
        {
            Phi = dlist[i].j0[edPhi];
            Psi = dlist[i].j0[edPsi];
            ca  = cb = co = ha = 0;
            for (j = 0; (j < nf); j++)
            {
                phi = dih[Phi][j];
                psi = dih[Psi][j];

                ca += interpolate(phi, psi, ca_sd);
                cb += interpolate(phi, psi, cb_sd);
                co += interpolate(phi, psi, co_sd);
                ha += interpolate(phi, psi, ha_sd);
            }
            fprintf(fp, "%12s  %10g  %10g  %10g  %10g\n",
                    dlist[i].name, ca/nf, ha/nf, co/nf, cb/nf);
        }
    }
    fprintf(fp, "\n");

    /* Free memory */
    done_shifts(ca_sd);
    done_shifts(cb_sd);
    done_shifts(co_sd);
    done_shifts(ha_sd);
}

int gmx_dos(int argc, char *argv[])
{
    const char         *desc[] = {
        "[TT]g_dos[tt] computes the Density of States from a simulations.",
        "In order for this to be meaningful the velocities must be saved",
        "in the trajecotry with sufficiently high frequency such as to cover",
        "all vibrations. For flexible systems that would be around a few fs",
        "between saving. Properties based on the DoS are printed on the",
        "standard output."
    };
    const char         *bugs[] = {
        "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)."
    };
    FILE               *fp, *fplog;
    t_topology          top;
    int                 ePBC = -1;
    t_trxframe          fr;
    matrix              box;
    int                 gnx;
    char                title[256];
    real                t0, t1, m;
    t_trxstatus        *status;
    int                 nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
    double              rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
    real              **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
    output_env_t        oenv;
    gmx_fft_t           fft;
    double              cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
    double              wCdiff, wSdiff, wAdiff, wEdiff;

    static     gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
    static     gmx_bool bRecip   = FALSE, bDump = FALSE;
    static     real     Temp     = 298.15, toler = 1e-6;
    t_pargs             pa[]     = {
        { "-v", FALSE, etBOOL, {&bVerbose},
          "Be loud and noisy." },
        { "-recip", FALSE, etBOOL, {&bRecip},
          "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
        { "-abs", FALSE, etBOOL, {&bAbsolute},
          "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
        { "-normdos", FALSE, etBOOL, {&bNormalize},
          "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
        { "-T", FALSE, etREAL, {&Temp},
          "Temperature in the simulation" },
        { "-toler", FALSE, etREAL, {&toler},
          "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
        { "-dump", FALSE, etBOOL, {&bDump},
          "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
    };

    t_filenm            fnm[] = {
        { efTRN, "-f",    NULL,    ffREAD  },
        { efTPX, "-s",    NULL,    ffREAD  },
        { efNDX, NULL,    NULL,    ffOPTRD },
        { efXVG, "-vacf", "vacf",  ffWRITE },
        { efXVG, "-mvacf", "mvacf", ffWRITE },
        { efXVG, "-dos",  "dos",   ffWRITE },
        { efLOG, "-g",    "dos",   ffWRITE },
    };
#define NFILE asize(fnm)
    int                 npargs;
    t_pargs            *ppa;
    const char         *DoSlegend[] = {
        "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
    };

    npargs = asize(pa);
    ppa    = add_acf_pargs(&npargs, pa);
    parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
                      NFILE, fnm, npargs, ppa, asize(desc), desc,
                      asize(bugs), bugs, &oenv);

    beta = 1/(Temp*BOLTZ);
    if (bDump)
    {
        printf("Dumping reference figures. Thanks for your patience.\n");
        dump_fy(oenv, toler);
        dump_w(oenv, beta);
        exit(0);
    }

    fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
    fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
    please_cite(fplog, "Pascal2011a");
    please_cite(fplog, "Caleman2011b");

    read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
                  TRUE);
    V     = det(box);
    tmass = 0;
    for (i = 0; (i < top.atoms.nr); i++)
    {
        tmass += top.atoms.atom[i].m;
    }

    Natom = top.atoms.nr;
    Nmol  = top.mols.nr;
    gnx   = Natom*DIM;

    /* Correlation stuff */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if (MULTISIM(cr) && MASTER(cr))
      {
	 if (MULTIMASTER(cr))
	 {
	    /* Log file is not yet available, so if there's a
	     * problem we can only write to stderr. */
	    fpmulti = stderr;
	 }
	 else
	 {
	    fpmulti = NULL;
	 }
	 check_multi_int(fpmulti, cr->ms, sim_part, "simulation part", TRUE);
      }
   }
   else
   {
      bAppendFiles = FALSE;
   }

   if (!bAppendFiles)
   {
      sim_part_fn = sim_part;
   }

   if (bAddPart)
   {
      /* Rename all output files (except checkpoint files) */
      /* create new part name first (zero-filled) */
      sprintf(suffix, "%s%04d", part_suffix, sim_part_fn);

      add_suffix_to_output_names(fnm, NFILE, suffix);
      if (MULTIMASTER(cr))
      {
	 fprintf(stdout, "Checkpoint file is from part %d, new output files will be suffixed '%s'.\n", sim_part-1, suffix);
      }
   }

   Flags = opt2bSet("-rerun", NFILE, fnm) ? MD_RERUN : 0;
   Flags = Flags | (bSepPot       ? MD_SEPPOT       : 0);
   Flags = Flags | (bDDBondCheck  ? MD_DDBONDCHECK  : 0);
   Flags = Flags | (bDDBondComm   ? MD_DDBONDCOMM   : 0);
   Flags = Flags | (bTunePME      ? MD_TUNEPME      : 0);
   Flags = Flags | (bTestVerlet   ? MD_TESTVERLET   : 0);
   Flags = Flags | (bConfout      ? MD_CONFOUT      : 0);
   Flags = Flags | (bRerunVSite   ? MD_RERUN_VSITE  : 0);
   Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
   Flags = Flags | (bAppendFiles  ? MD_APPENDFILES  : 0);
   Flags = Flags | (opt2parg_bSet("-append", asize(pa), pa) ? MD_APPENDFILESSET : 0);
   Flags = Flags | (bKeepAndNumCPT ? MD_KEEPANDNUMCPT : 0);
   Flags = Flags | (sim_part > 1    ? MD_STARTFROMCPT : 0);
   Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
   Flags = Flags | (bIMDwait      ? MD_IMDWAIT      : 0);
   Flags = Flags | (bIMDterm      ? MD_IMDTERM      : 0);
   Flags = Flags | (bIMDpull      ? MD_IMDPULL      : 0);
   Flags = Flags | ((opt2bSet("-at", NFILE, fnm)) ? MD_ADAPTIVETEMPERING : 0);
   Flags = Flags | ((opt2bSet("-addtop", NFILE, fnm)) ? MD_MULTOP : 0);

   /* We postpone opening the log file if we are appending, so we can
      first truncate the old log file and append to the correct position
      there instead.  */
   if ((MASTER(cr) || bSepPot) && !bAppendFiles)
   {
      gmx_log_open(ftp2fn(efLOG, NFILE, fnm), cr,
	    !bSepPot, Flags & MD_APPENDFILES, &fplog);
      please_cite(fplog, "Hess2008b");
      please_cite(fplog, "Spoel2005a");
      please_cite(fplog, "Lindahl2001a");
      please_cite(fplog, "Berendsen95a");
   }
   else if (!MASTER(cr) && bSepPot)
   {
      gmx_log_open(ftp2fn(efLOG, NFILE, fnm), cr, !bSepPot, Flags, &fplog);
   }
   else
   {
      fplog = NULL;
   }

   ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
   ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
   ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);

   rc = mdrunner(&hw_opt, fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
	 nstglobalcomm, ddxyz, dd_node_order, rdd, rconstr,
	 dddlb_opt[0], dlb_scale, ddcsx, ddcsy, ddcsz,
	 nbpu_opt[0], nstlist,
	 nsteps, nstepout, resetstep,
	 nmultisim, repl_ex_nst, repl_ex_nex, repl_ex_seed,
	 pforce, cpt_period, max_hours, deviceOptions, imdport, Flags);

   /* Log file has to be closed in mdrunner if we are appending to it
      (fplog not set here) */
   if (MASTER(cr) && !bAppendFiles)
   {
      gmx_log_close(fplog);
   }

   return rc;
}

//.........这里部分代码省略.........
        }
        pr->grn      = prframecurrent->grn;
        pr->binwidth = prframecurrent->binwidth;
        /* summ up gr and fill r */
        for (i = 0; i < prframecurrent->grn; i++)
        {
            pr->gr[i] += prframecurrent->gr[i];
            pr->r[i]   = prframecurrent->r[i];
        }
        /* normalize histo */
        normalize_probability(prframecurrent->grn, prframecurrent->gr);
        /* convert p(r) to sq */
        sqframecurrent = convert_histogram_to_intensity_curve(prframecurrent, start_q, end_q, q_step);
        /* print frame data if needed */
        if (opt2fn_null("-prframe", NFILE, fnm))
        {
            snew(hdr, 25);
            snew(suffix, GMX_PATH_MAX);
            /* prepare header */
            sprintf(hdr, "g(r), t = %f", t);
            /* prepare output filename */
            fnmdup = dup_tfn(NFILE, fnm);
            sprintf(suffix, "-t%.2f", t);
            add_suffix_to_output_names(fnmdup, NFILE, suffix);
            fp = xvgropen(opt2fn_null("-prframe", NFILE, fnmdup), hdr, "Distance (nm)", "Probability", oenv);
            for (i = 0; i < prframecurrent->grn; i++)
            {
                fprintf(fp, "%10.6f%10.6f\n", prframecurrent->r[i], prframecurrent->gr[i]);
            }
            done_filenms(NFILE, fnmdup);
            fclose(fp);
            sfree(hdr);
            sfree(suffix);
            sfree(fnmdup);
        }
        if (opt2fn_null("-sqframe", NFILE, fnm))
        {
            snew(hdr, 25);
            snew(suffix, GMX_PATH_MAX);
            /* prepare header */
            sprintf(hdr, "I(q), t = %f", t);
            /* prepare output filename */
            fnmdup = dup_tfn(NFILE, fnm);
            sprintf(suffix, "-t%.2f", t);
            add_suffix_to_output_names(fnmdup, NFILE, suffix);
            fp = xvgropen(opt2fn_null("-sqframe", NFILE, fnmdup), hdr, "q (nm^-1)", "s(q)/s(0)", oenv);
            for (i = 0; i < sqframecurrent->qn; i++)
            {
                fprintf(fp, "%10.6f%10.6f\n", sqframecurrent->q[i], sqframecurrent->s[i]);
            }
            done_filenms(NFILE, fnmdup);
            fclose(fp);
            sfree(hdr);
            sfree(suffix);
            sfree(fnmdup);
        }
        /* free pr structure */
        sfree(prframecurrent->gr);
        sfree(prframecurrent->r);
        sfree(prframecurrent);
        /* free sq structure */
        sfree(sqframecurrent->q);
        sfree(sqframecurrent->s);
        sfree(sqframecurrent);
    }
    while (read_next_x(oenv, status, &t, x, box));
    close_trj(status);

    /* normalize histo */
    normalize_probability(pr->grn, pr->gr);
    sq = convert_histogram_to_intensity_curve(pr, start_q, end_q, q_step);
    /* prepare pr.xvg */
    fp = xvgropen(opt2fn_null("-pr", NFILE, fnm), "G(r)", "Distance (nm)", "Probability", oenv);
    for (i = 0; i < pr->grn; i++)
    {
        fprintf(fp, "%10.6f%10.6f\n", pr->r[i], pr->gr[i]);
    }
    xvgrclose(fp);

    /* prepare sq.xvg */
    fp = xvgropen(opt2fn_null("-sq", NFILE, fnm), "I(q)", "q (nm^-1)", "s(q)/s(0)", oenv);
    for (i = 0; i < sq->qn; i++)
    {
        fprintf(fp, "%10.6f%10.6f\n", sq->q[i], sq->s[i]);
    }
    xvgrclose(fp);
    /*
     * Clean up memory
     */
    sfree(pr->gr);
    sfree(pr->r);
    sfree(pr);
    sfree(sq->q);
    sfree(sq->s);
    sfree(sq);

    please_cite(stdout, "Garmay2012");

    return 0;
}

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

    npargs = asize(pa);
    ppa    = add_acf_pargs(&npargs, pa);
    parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
                      NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs,
                      &oenv);

    mult   = 4;
    maxang = 360.0;
    bRb    = FALSE;
    switch (opt[0][0])
    {
    case 'a':
        mult   = 3;
        maxang = 180.0;
        break;
    case 'd':
        break;
    case 'i':
        break;
    case 'r':
        bRb = TRUE;
        break;
    }

    if (opt2bSet("-or", NFILE, fnm))
    {
        if (mult != 4)
        {
            gmx_fatal(FARGS, "Can not combine angles with trn dump");
        }
        else
        {
            please_cite(stdout, "Mu2005a");
        }
    }

    /* Calculate bin size */
    maxangstat = (int)(maxang/binwidth+0.5);
    binwidth   = maxang/maxangstat;

    rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
    nangles = isize/mult;
    if ((isize % mult) != 0)
    {
        gmx_fatal(FARGS, "number of index elements not multiple of %d, "
                  "these can not be %s\n",
                  mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
    }


    /* Check whether specific analysis has to be performed */
    bCorr  = opt2bSet("-oc", NFILE, fnm);
    bAver  = opt2bSet("-ov", NFILE, fnm);
    bTrans = opt2bSet("-ot", NFILE, fnm);
    bFrac  = opt2bSet("-of", NFILE, fnm);
    if (bTrans && opt[0][0] != 'd')
    {
        fprintf(stderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
        bTrans = FALSE;
    }

    if (bChandler && !bCorr)
    {
        bCorr = TRUE;
    }

//.........这里部分代码省略.........
    atom_id     *ind_fit, **ind_rms, *ind_m = NULL, *rev_ind_m = NULL, *ind_rms_m =
        NULL;
    char        *gn_fit, **gn_rms;
    t_rgb        rlo, rhi;
    output_env_t oenv;
    gmx_rmpbc_t  gpbc = NULL;

    t_filenm     fnm[] =
    {
        { efTPS, NULL, NULL, ffREAD },
        { efTRX, "-f", NULL, ffREAD },
        { efTRX, "-f2", NULL, ffOPTRD },
        { efNDX, NULL, NULL, ffOPTRD },
        { efXVG, NULL, "rmsd", ffWRITE },
        { efXVG, "-mir", "rmsdmir", ffOPTWR },
        { efXVG, "-a", "avgrp", ffOPTWR },
        { efXVG, "-dist", "rmsd-dist", ffOPTWR },
        { efXPM, "-m", "rmsd", ffOPTWR },
        { efDAT, "-bin", "rmsd", ffOPTWR },
        { efXPM, "-bm", "bond", ffOPTWR }
    };
#define NFILE asize(fnm)

    if (!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))
    {
        return 0;
    }
    /* parse enumerated options: */
    ewhat = nenum(what);
    if (ewhat == ewRho || ewhat == ewRhoSc)
    {
        please_cite(stdout, "Maiorov95");
    }
    efit   = nenum(fit);
    bFit   = efit == efFit;
    bReset = efit == efReset;
    if (bFit)
    {
        bReset = TRUE; /* for fit, reset *must* be set */
    }
    else
    {
        bFitAll = FALSE;
    }

    /* mark active cmdline options */
    bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
    bFile2  = opt2bSet("-f2", NFILE, fnm);
    bMat    = opt2bSet("-m", NFILE, fnm);
    bBond   = opt2bSet("-bm", NFILE, fnm);
    bDelta  = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
                                 *	your RMSD matrix (hidden option       */
    bNorm   = opt2bSet("-a", NFILE, fnm);
    bFreq2  = opt2parg_bSet("-skip2", asize(pa), pa);
    if (freq <= 0)
    {
        fprintf(stderr, "The number of frames to skip is <= 0. "
                "Writing out all frames.\n\n");
        freq = 1;
    }
    if (!bFreq2)
    {
        freq2 = freq;
    }

void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
               gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
               gmx_bool bOutFile, unsigned long Flags)
{
    t_pull       *pull;
    t_pull_group *pgrp;
    int           c, g, start = 0, end = 0, m;

    pull = ir->pull;

    pull->ePBC = ir->ePBC;
    switch (pull->ePBC)
    {
        case epbcNONE: pull->npbcdim = 0; break;
        case epbcXY:   pull->npbcdim = 2; break;
        default:       pull->npbcdim = 3; break;
    }

    if (fplog)
    {
        gmx_bool bAbs, bCos;

        bAbs = FALSE;
        for (c = 0; c < pull->ncoord; c++)
        {
            if (pull->group[pull->coord[c].group[0]].nat == 0 ||
                pull->group[pull->coord[c].group[1]].nat == 0)
            {
                bAbs = TRUE;
            }
        }
        
        fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
                EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
        fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
                pull->ncoord, pull->ncoord == 1 ? "" : "s",
                pull->ngroup, pull->ngroup == 1 ? "" : "s");
        if (bAbs)
        {
            fprintf(fplog, "with an absolute reference\n");
        }
        bCos = FALSE;
        for (g = 0; g < pull->ngroup; g++)
        {
            if (pull->group[g].nat > 1 &&
                pull->group[g].pbcatom < 0)
            {
                /* We are using cosine weighting */
                fprintf(fplog, "Cosine weighting is used for group %d\n", g);
                bCos = TRUE;
            }
        }
        if (bCos)
        {
            please_cite(fplog, "Engin2010");
        }
    }

    /* We always add the virial contribution,
     * except for geometry = direction_periodic where this is impossible.
     */
    pull->bVirial = (pull->eGeom != epullgDIRPBC);
    if (getenv("GMX_NO_PULLVIR") != NULL)
    {
        if (fplog)
        {
            fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
        }
        pull->bVirial = FALSE;
    }

    if (cr && PARTDECOMP(cr))
    {
        pd_at_range(cr, &start, &end);
    }
    pull->rbuf     = NULL;
    pull->dbuf     = NULL;
    pull->dbuf_cyl = NULL;
    pull->bRefAt   = FALSE;
    pull->cosdim   = -1;
    for (g = 0; g < pull->ngroup; g++)
    {
        pgrp           = &pull->group[g];
        pgrp->epgrppbc = epgrppbcNONE;
        if (pgrp->nat > 0)
        {
            /* Determine if we need to take PBC into account for calculating
             * the COM's of the pull groups.
             */
            for (m = 0; m < pull->npbcdim; m++)
            {
                if (pull->dim[m] && pgrp->nat > 1)
                {
                    if (pgrp->pbcatom >= 0)
                    {
                        pgrp->epgrppbc = epgrppbcREFAT;
                        pull->bRefAt   = TRUE;
                    }
                    else
                    {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            filename = ftp2fn(efTRN, nfile, fnm);
            switch (fn2ftp(filename))
            {
                case efTRR:
                case efTRN:
                    of->fp_trn = gmx_trr_open(filename, filemode);
                    break;
                case efTNG:
                    gmx_tng_open(filename, filemode[0], &of->tng);
                    if (filemode[0] == 'w')
                    {
                        gmx_tng_prepare_md_writing(of->tng, top_global, ir);
                    }
                    bCiteTng = TRUE;
                    break;
                default:
                    gmx_incons("Invalid full precision file format");
            }
        }
        if (EI_DYNAMICS(ir->eI) &&
            ir->nstxout_compressed > 0)
        {
            const char *filename;
            filename = ftp2fn(efCOMPRESSED, nfile, fnm);
            switch (fn2ftp(filename))
            {
                case efXTC:
                    of->fp_xtc                  = open_xtc(filename, filemode);
                    break;
                case efTNG:
                    gmx_tng_open(filename, filemode[0], &of->tng_low_prec);
                    if (filemode[0] == 'w')
                    {
                        gmx_tng_prepare_low_prec_writing(of->tng_low_prec, top_global, ir);
                    }
                    bCiteTng = TRUE;
                    break;
                default:
                    gmx_incons("Invalid reduced precision file format");
            }
        }
        if (EI_DYNAMICS(ir->eI) || EI_ENERGY_MINIMIZATION(ir->eI))
        {
            of->fp_ene = open_enx(ftp2fn(efEDR, nfile, fnm), filemode);
        }
        of->fn_cpt = opt2fn("-cpo", nfile, fnm);

        if ((ir->efep != efepNO || ir->bSimTemp) && ir->fepvals->nstdhdl > 0 &&
            (ir->fepvals->separate_dhdl_file == esepdhdlfileYES ) &&
            EI_DYNAMICS(ir->eI))
        {
            if (bAppendFiles)
            {
                of->fp_dhdl = gmx_fio_fopen(opt2fn("-dhdl", nfile, fnm), filemode);
            }
            else
            {
                of->fp_dhdl = open_dhdl(opt2fn("-dhdl", nfile, fnm), ir, oenv);
            }
        }

        if (opt2bSet("-field", nfile, fnm) &&
            (ir->ex[XX].n || ir->ex[YY].n || ir->ex[ZZ].n))
        {
            if (bAppendFiles)
            {
                of->fp_field = gmx_fio_fopen(opt2fn("-field", nfile, fnm),
                                             filemode);
            }
            else
            {
                of->fp_field = xvgropen(opt2fn("-field", nfile, fnm),
                                        "Applied electric field", "Time (ps)",
                                        "E (V/nm)", oenv);
            }
        }

        /* Set up atom counts so they can be passed to actual
           trajectory-writing routines later. Also, XTC writing needs
           to know what (and how many) atoms might be in the XTC
           groups, and how to look up later which ones they are. */
        of->natoms_global       = top_global->natoms;
        of->groups              = &top_global->groups;
        of->natoms_x_compressed = 0;
        for (i = 0; (i < top_global->natoms); i++)
        {
            if (ggrpnr(of->groups, egcCompressedX, i) == 0)
            {
                of->natoms_x_compressed++;
            }
        }
    }

    if (bCiteTng)
    {
        please_cite(fplog, "Lundborg2014");
    }

    return of;
}

void init_md(FILE *fplog,
			 t_commrec *cr,t_inputrec *ir,real *t,real *t0,
			 real *lambda,real *lam0,
			 t_nrnb *nrnb,gmx_mtop_t *mtop,
			 gmx_stochd_t *sd,
			 int nfile,t_filenm fnm[],
			 int *fp_trn,int *fp_xtc,int *fp_ene,char **fn_cpt,
			 FILE **fp_dgdl,FILE **fp_field,
			 t_mdebin **mdebin,
			 tensor force_vir,tensor shake_vir,rvec mu_tot,
			 bool *bNEMD,bool *bSimAnn,t_vcm **vcm, unsigned long Flags)
{
  int  i,j,n;
  real tmpt,mod;
  char filemode[2];

  sprintf(filemode, (Flags & MD_APPENDFILES) ? "a" : "w");
	
  /* Initial values */
  *t = *t0       = ir->init_t;
  if (ir->efep != efepNO) {
    *lam0 = ir->init_lambda;
    *lambda = *lam0 + ir->init_step*ir->delta_lambda;
  }
  else {
    *lambda = *lam0   = 0.0;
  } 
 
  *bSimAnn=FALSE;
  for(i=0;i<ir->opts.ngtc;i++) {
    /* set bSimAnn if any group is being annealed */
    if(ir->opts.annealing[i]!=eannNO)
      *bSimAnn = TRUE;
  }
  if (*bSimAnn) {
    update_annealing_target_temp(&(ir->opts),ir->init_t);
  }

  *bNEMD = (ir->opts.ngacc > 1) || (norm(ir->opts.acc[0]) > 0);
  
  if (sd && (ir->eI == eiBD || EI_SD(ir->eI) || ir->etc == etcVRESCALE)) {
    *sd = init_stochd(fplog,ir);
  }

  if (vcm) {
    *vcm = init_vcm(fplog,&mtop->groups,ir);
  }
   
  if (EI_DYNAMICS(ir->eI) && !(Flags & MD_APPENDFILES)) {
    if (ir->etc == etcBERENDSEN) {
      please_cite(fplog,"Berendsen84a");
    }
    if (ir->etc == etcVRESCALE) {
      please_cite(fplog,"Bussi2007a");
    }
  }
 
  init_nrnb(nrnb);

  if (nfile != -1)
  {
	  *fp_trn = -1;
	  *fp_ene = -1;
	  *fp_xtc = -1;
	  
	  if (MASTER(cr)) 
	  {
		  *fp_trn = open_trn(ftp2fn(efTRN,nfile,fnm), filemode);
		  if (ir->nstxtcout > 0)
		  {
			  *fp_xtc = open_xtc(ftp2fn(efXTC,nfile,fnm), filemode);
		  }
		  *fp_ene = open_enx(ftp2fn(efENX,nfile,fnm), filemode);
		  *fn_cpt = opt2fn("-cpo",nfile,fnm);
		  
		  if ((fp_dgdl != NULL) && ir->efep!=efepNO)
		  {
			  if(Flags & MD_APPENDFILES)
			  {
				  *fp_dgdl= gmx_fio_fopen(opt2fn("-dgdl",nfile,fnm),filemode);
			  }
			  else
			  {
				  *fp_dgdl = xvgropen(opt2fn("-dgdl",nfile,fnm),
									  "dG/d\\8l\\4","Time (ps)",
									  "dG/d\\8l\\4 (kJ mol\\S-1\\N [\\8l\\4]\\S-1\\N)");
			  }
		  }
		  
		  if ((fp_field != NULL) && (ir->ex[XX].n || ir->ex[YY].n ||ir->ex[ZZ].n))
		  {
			  if(Flags & MD_APPENDFILES)
			  {
				  *fp_field = gmx_fio_fopen(opt2fn("-field",nfile,fnm),filemode);
			  }
			  else
			  {				  
				  *fp_field = xvgropen(opt2fn("-field",nfile,fnm),
									   "Applied electric field","Time (ps)",
									   "E (V/nm)");
//.........这里部分代码省略.........

int gmx_dielectric(int argc, char *argv[])
{
    const char  *desc[] = {
        "[THISMODULE] calculates frequency dependent dielectric constants",
        "from the autocorrelation function of the total dipole moment in",
        "your simulation. This ACF can be generated by [gmx-dipoles].",
        "The functional forms of the available functions are:[PAR]",
        "One parameter:    y = [EXP]-a[SUB]1[sub] x[exp],[BR]",
        "Two parameters:   y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp],[BR]",
        "Three parameters: y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp] + (1 - a[SUB]2[sub]) [EXP]-a[SUB]3[sub] x[exp].[BR]",
        "Start values for the fit procedure can be given on the command line.",
        "It is also possible to fix parameters at their start value, use [TT]-fix[tt]",
        "with the number of the parameter you want to fix.",
        "[PAR]",
        "Three output files are generated, the first contains the ACF,",
        "an exponential fit to it with 1, 2 or 3 parameters, and the",
        "numerical derivative of the combination data/fit.",
        "The second file contains the real and imaginary parts of the",
        "frequency-dependent dielectric constant, the last gives a plot",
        "known as the Cole-Cole plot, in which the imaginary",
        "component is plotted as a function of the real component.",
        "For a pure exponential relaxation (Debye relaxation) the latter",
        "plot should be one half of a circle."
    };
    t_filenm     fnm[] = {
        { efXVG, "-f", "dipcorr", ffREAD  },
        { efXVG, "-d", "deriv",  ffWRITE },
        { efXVG, "-o", "epsw",   ffWRITE },
        { efXVG, "-c", "cole",   ffWRITE }
    };
#define NFILE asize(fnm)
    output_env_t oenv;
    int          i, j, nx, ny, nxtail, eFitFn, nfitparm;
    real         dt, integral, fitintegral, *fitparms, fac, rffac;
    double     **yd;
    real       **y;
    const char  *legend[] = { "Correlation", "Std. Dev.", "Fit", "Combined", "Derivative" };
    static int   fix      = 0, bFour = 0, bX = 1, nsmooth = 3;
    static real  tendInt  = 5.0, tbegin = 5.0, tend = 500.0;
    static real  A        = 0.5, tau1 = 10.0, tau2 = 1.0, eps0 = 80, epsRF = 78.5, tail = 500.0;
    real         lambda;
    t_pargs      pa[] = {
        { "-fft", FALSE, etBOOL, {&bFour},
          "use fast fourier transform for correlation function" },
        { "-x1",  FALSE, etBOOL, {&bX},
          "use first column as [IT]x[it]-axis rather than first data set" },
        { "-eint", FALSE, etREAL, {&tendInt},
          "Time to end the integration of the data and start to use the fit"},
        { "-bfit", FALSE, etREAL, {&tbegin},
          "Begin time of fit" },
        { "-efit", FALSE, etREAL, {&tend},
          "End time of fit" },
        { "-tail", FALSE, etREAL, {&tail},
          "Length of function including data and tail from fit" },
        { "-A", FALSE, etREAL, {&A},
          "Start value for fit parameter A" },
        { "-tau1", FALSE, etREAL, {&tau1},
          "Start value for fit parameter [GRK]tau[grk]1" },
        { "-tau2", FALSE, etREAL, {&tau2},
          "Start value for fit parameter [GRK]tau[grk]2" },
        { "-eps0", FALSE, etREAL, {&eps0},
          "[GRK]epsilon[grk]0 of your liquid" },
        { "-epsRF", FALSE, etREAL, {&epsRF},
          "[GRK]epsilon[grk] of the reaction field used in your simulation. A value of 0 means infinity." },
        { "-fix", FALSE, etINT,  {&fix},
          "Fix parameters at their start values, A (2), tau1 (1), or tau2 (4)" },
        { "-ffn",    FALSE, etENUM, {s_ffn},
          "Fit function" },
        { "-nsmooth", FALSE, etINT, {&nsmooth},
          "Number of points for smoothing" }
    };

    if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
    {
        return 0;
    }
    please_cite(stdout, "Spoel98a");
    printf("WARNING: non-polarizable models can never yield an infinite\n"
           "dielectric constant that is different from 1. This is incorrect\n"
           "in the reference given above (Spoel98a).\n\n");


    nx     = read_xvg(opt2fn("-f", NFILE, fnm), &yd, &ny);
    dt     = yd[0][1] - yd[0][0];
    nxtail = min(tail/dt, nx);

    printf("Read data set containing %d colums and %d rows\n", ny, nx);
    printf("Assuming (from data) that timestep is %g, nxtail = %d\n",
           dt, nxtail);
    snew(y, 6);
    for (i = 0; (i < ny); i++)
    {
        snew(y[i], max(nx, nxtail));
    }
    for (i = 0; (i < nx); i++)
    {
        y[0][i] = yd[0][i];
        for (j = 1; (j < ny); j++)
        {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        { "-sum",     FALSE, etBOOL, {&bSum},
          "Average folding before computing [GRK]chi[grk]^2" },
        { "-discrete", FALSE, etBOOL, {&bDiscrete},
          "Use a discrete folding criterion (F <-> U) or a continuous one" },
        { "-mult",    FALSE, etINT, {&nmult},
          "Factor to multiply the data with before discretization" }
    };
#define NPA asize(pa)

    FILE        *fp;
    real         dt_t, dt_d, dt_d2;
    int          nset_t, nset_d, nset_d2, n_t, n_d, n_d2, i;
    const char  *tfile, *dfile, *dfile2;
    t_remd_data  remd;
    output_env_t oenv;

    t_filenm     fnm[] = {
        { efXVG, "-f",    "temp",    ffREAD   },
        { efXVG, "-d",    "data",    ffREAD   },
        { efXVG, "-d2",   "data2",   ffOPTRD  },
        { efXVG, "-o",    "ft_all",  ffWRITE  },
        { efXVG, "-o2",   "it_all",  ffOPTWR  },
        { efXVG, "-o3",   "ft_repl", ffOPTWR  },
        { efXVG, "-ee",   "err_est", ffOPTWR  },
        { efLOG, "-g",    "remd",    ffWRITE  },
        { efXVG, "-m",    "melt",    ffWRITE  }
    };
#define NFILE asize(fnm)

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

#ifdef HAVE_LIBGSL
    please_cite(stdout, "Spoel2006d");
    if (cutoff < 0)
    {
        gmx_fatal(FARGS, "cutoff should be >= 0 (rather than %f)", cutoff);
    }

    tfile   = opt2fn("-f", NFILE, fnm);
    dfile   = opt2fn("-d", NFILE, fnm);
    dfile2  = opt2fn_null("-d2", NFILE, fnm);

    fp = ffopen(opt2fn("-g", NFILE, fnm), "w");

    remd.temp = read_xvg_time(tfile, bHaveT,
                              opt2parg_bSet("-b", NPA, pa), tb,
                              opt2parg_bSet("-e", NPA, pa), te,
                              nreplica, &nset_t, &n_t, &dt_t, &remd.time);
    printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_t, n_t, tfile, dt_t);
    sfree(remd.time);

    remd.data = read_xvg_time(dfile, bHaveT,
                              opt2parg_bSet("-b", NPA, pa), tb,
                              opt2parg_bSet("-e", NPA, pa), te,
                              nreplica, &nset_d, &n_d, &dt_d, &remd.time);
    printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_d, n_d, dfile, dt_d);

    if ((nset_t != nset_d) || (n_t != n_d) || (dt_t != dt_d))
    {
        gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file",
                  tfile, dfile);
    }

    if (dfile2)
    {