static void dump_slab_dipoles(const char *fn,int idim,int nslice,
                              rvec slab_dipole[], matrix box,int nframes,
                              const output_env_t oenv)
{
    FILE *fp;
    char buf[STRLEN];
    int  i;
    real mutot;
    const char *leg_dim[4] = { 
        "\\f{12}m\\f{4}\\sX\\N",
        "\\f{12}m\\f{4}\\sY\\N",
        "\\f{12}m\\f{4}\\sZ\\N",
        "\\f{12}m\\f{4}\\stot\\N"
    };
  
    sprintf(buf,"Box-%c (nm)",'X'+idim);
    fp = xvgropen(fn,"Average dipole moment per slab",buf,"\\f{12}m\\f{4} (D)",
                  oenv);
    xvgr_legend(fp,DIM,leg_dim,oenv); 
    for(i=0; (i<nslice); i++) {
        mutot = norm(slab_dipole[i])/nframes;
        fprintf(fp,"%10.3f  %10.3f  %10.3f  %10.3f  %10.3f\n",
                ((i+0.5)*box[idim][idim])/nslice,
                slab_dipole[i][XX]/nframes,
                slab_dipole[i][YY]/nframes,
                slab_dipole[i][ZZ]/nframes,
                mutot);
    }
    ffclose(fp);
    do_view(oenv,fn,"-autoscale xy -nxy");
}

void view_all(const output_env_t oenv,int nf, t_filenm fnm[])
{
  int i;
  
  for(i=0; i<nf; i++)
    if ( can_view(fnm[i].ftp) && is_output(&(fnm[i])) && 
	 ( ! is_optional(&(fnm[i])) || is_set(&(fnm[i])) ) )
      do_view(oenv,fnm[i].fns[0], NULL);
}

void view_all(int nf, t_filenm fnm[])
{
  int i;
  
  for(i=0; i<nf; i++)
    if ( can_view(fnm[i].ftp) && is_output(&(fnm[i])) && 
	 ( ! is_optional(&(fnm[i])) || is_set(&(fnm[i])) ) )
      do_view(fnm[i].fn, NULL);
}

void view_all(const gmx_output_env_t *oenv, int nf, t_filenm fnm[])
{
    int i;

    for (i = 0; i < nf; i++)
    {
        if (can_view(fnm[i].ftp) && is_output(&(fnm[i])) &&
            ( !is_optional(&(fnm[i])) || is_set(&(fnm[i])) ) )
        {
            do_view(oenv, fnm[i].filenames[0].c_str(), nullptr);
        }
    }
}

int gmx_rama(int argc,char *argv[])
{
  const char *desc[] = {
    "[TT]g_rama[tt] selects the [GRK]phi[grk]/[GRK]psi[grk] dihedral combinations from your topology file",
    "and computes these as a function of time.",
    "Using simple Unix tools such as [IT]grep[it] you can select out", 
    "specific residues."
  };

  FILE      *out;
  t_xrama   *xr;
  int       j;
  output_env_t oenv;
  t_filenm  fnm[] = {
    { efTRX, "-f", NULL,  ffREAD },
    { efTPX, NULL, NULL,  ffREAD },
    { efXVG, NULL, "rama",ffWRITE }
  };
#define NFILE asize(fnm)

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

		      
  snew(xr,1);
  init_rama(oenv,ftp2fn(efTRX,NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),xr,3);
  
  out=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Ramachandran Plot","Phi","Psi",oenv);
  xvgr_line_props(out,0,elNone,ecFrank,oenv);
  xvgr_view(out,0.2,0.2,0.8,0.8,oenv);
  xvgr_world(out,-180,-180,180,180,oenv);
  fprintf(out,"@    xaxis  tick on\[email protected]    xaxis  tick major 60\[email protected]    xaxis  tick minor 30\n");
  fprintf(out,"@    yaxis  tick on\[email protected]    yaxis  tick major 60\[email protected]    yaxis  tick minor 30\n");
  fprintf(out,"@ s0 symbol 2\[email protected] s0 symbol size 0.4\[email protected] s0 symbol fill 1\n");
  
  j=0;
  do {
    plot_rama(out,xr);
    j++;
  } while (new_data(xr));
  fprintf(stderr,"\n");
  ffclose(out);
  
  do_view(oenv,ftp2fn(efXVG,NFILE,fnm),NULL);
  
  thanx(stderr);
  
  return 0;
}

//.........这里部分代码省略.........
                  whatxvglabel[ewhat], oenv);
    if (output_env_get_print_xvgr_codes(oenv))
    {
        fprintf(fp, "@ subtitle \"%s%s after %s%s%s\"\n",
                (nrms == 1) ? "" : "of ", gn_rms[0], fitgraphlabel[efit],
                bFit     ? " to " : "", bFit ? gn_fit : "");
    }
    if (nrms != 1)
    {
        xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
    }
    for (i = 0; (i < teller); i++)
    {
        if (bSplit && i > 0 &&
            abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
        {
            fprintf(fp, "&\n");
        }
        fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]);
        for (j = 0; (j < nrms); j++)
        {
            fprintf(fp, " %12.7f", rls[j][i]);
            if (bAv)
            {
                rlstot += rls[j][i];
            }
        }
        fprintf(fp, "\n");
    }
    gmx_ffclose(fp);

    if (bMirror)
    {
        /* Write the mirror RMSD's to file */
        sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);
        sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);
        fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
                      buf2, oenv);
        if (nrms == 1)
        {
            if (output_env_get_print_xvgr_codes(oenv))
            {
                fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
                        gn_rms[0], gn_fit);
            }
        }
        else
        {
            if (output_env_get_print_xvgr_codes(oenv))
            {
                fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", gn_fit);
            }
            xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
        }
        for (i = 0; (i < teller); i++)
        {
            if (bSplit && i > 0 && abs(time[i]) < 1e-5)
            {
                fprintf(fp, "&\n");
            }
            fprintf(fp, "%12.7f", time[i]);
            for (j = 0; (j < nrms); j++)
            {
                fprintf(fp, " %12.7f", rlsm[j][i]);
            }
            fprintf(fp, "\n");
        }
        gmx_ffclose(fp);
    }

    if (bAv)
    {
        sprintf(buf, "Average %s", whatxvgname[ewhat]);
        sprintf(buf2, "Average %s", whatxvglabel[ewhat]);
        fp = xvgropen(opt2fn("-a", NFILE, fnm), buf, "Residue", buf2, oenv);
        for (j = 0; (j < nrms); j++)
        {
            fprintf(fp, "%10d  %10g\n", j, rlstot/teller);
        }
        gmx_ffclose(fp);
    }

    if (bNorm)
    {
        fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv);
        for (j = 0; (j < irms[0]); j++)
        {
            fprintf(fp, "%10d  %10g\n", j, rlsnorm[j]/teller);
        }
        gmx_ffclose(fp);
    }
    do_view(oenv, opt2fn_null("-a", NFILE, fnm), "-graphtype bar");
    do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
    do_view(oenv, opt2fn_null("-mir", NFILE, fnm), NULL);
    do_view(oenv, opt2fn_null("-m", NFILE, fnm), NULL);
    do_view(oenv, opt2fn_null("-bm", NFILE, fnm), NULL);
    do_view(oenv, opt2fn_null("-dist", NFILE, fnm), NULL);

    return 0;
}

//.........这里部分代码省略.........
            }
            if (bChandler)
            {
                mode = eacNormal;
            }
            else
            {
                mode = eacCos;
            }
            do_autocorr(opt2fn("-oc", NFILE, fnm), oenv,
                        "Dihedral Autocorrelation Function",
                        nframes, nangles, dih, dt, mode, bAverCorr);
        }
    }


    /* Determine the non-zero part of the distribution */
    for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
    {
        ;
    }
    for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
    {
        ;
    }

    aver = aver2 = 0;
    for (i = 0; (i < nframes); i++)
    {
        aver  += RAD2DEG*aver_angle[i];
        aver2 += sqr(RAD2DEG*aver_angle[i]);
    }
    aver   /= (real) nframes;
    aver2  /= (real) nframes;
    aversig = sqrt(aver2-sqr(aver));
    printf("Found points in the range from %d to %d (max %d)\n",
           first, last, maxangstat);
    printf(" < angle >  = %g\n", aver);
    printf("< angle^2 > = %g\n", aver2);
    printf("Std. Dev.   = %g\n", aversig);

    if (mult == 3)
    {
        sprintf(title, "Angle Distribution: %s", grpname);
    }
    else
    {
        sprintf(title, "Dihedral Distribution: %s", grpname);

        calc_distribution_props(maxangstat, angstat, -180.0, 0, NULL, &S2);
        fprintf(stderr, "Order parameter S^2 = %g\n", S2);
    }

    bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);

    out = xvgropen(opt2fn("-od", NFILE, fnm), title, "Degrees", "", oenv);
    if (output_env_get_print_xvgr_codes(oenv))
    {
        fprintf(out, "@    subtitle \"average angle: %g\\So\\N\"\n", aver);
    }
    norm_fac = 1.0/(nangles*nframes*binwidth);
    if (bPeriodic)
    {
        maxstat = 0;
        for (i = first; (i <= last); i++)
        {
            maxstat = max(maxstat, angstat[i]*norm_fac);
        }
        fprintf(out, "@with g0\n");
        fprintf(out, "@    world xmin -180\n");
        fprintf(out, "@    world xmax  180\n");
        fprintf(out, "@    world ymin 0\n");
        fprintf(out, "@    world ymax %g\n", maxstat*1.05);
        fprintf(out, "@    xaxis  tick major 60\n");
        fprintf(out, "@    xaxis  tick minor 30\n");
        fprintf(out, "@    yaxis  tick major 0.005\n");
        fprintf(out, "@    yaxis  tick minor 0.0025\n");
    }
    for (i = first; (i <= last); i++)
    {
        fprintf(out, "%10g  %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
    }
    if (bPeriodic)
    {
        /* print first bin again as last one */
        fprintf(out, "%10g  %10f\n", 180.0, angstat[0]*norm_fac);
    }

    ffclose(out);

    do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
    if (bAver)
    {
        do_view(oenv, opt2fn("-ov", NFILE, fnm), "-nxy");
    }

    thanx(stderr);

    return 0;
}

//.........这里部分代码省略.........
            dos[DOS][j] = bfac*dos[DOS][2*j];
        }
    }
    /* Normalize it */
    dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
    if (bNormalizeDos)
    {
        for (j = 0; (j < nframes/4); j++)
        {
            dos[DOS][j] *= 3*Natom/dostot;
        }
    }

    /* Now analyze it */
    DoS0 = dos[DOS][0];

    /* Note this eqn. is incorrect in Pascal2011a! */
    Delta = ((2*DoS0/(9*Natom))*std::sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
             std::pow((Natom/V), 1.0/3.0)*std::pow(6.0/M_PI, 2.0/3.0));
    f     = calc_fluidicity(Delta, toler);
    y     = calc_y(f, Delta, toler);
    z     = calc_compress(y);
    Sig   = BOLTZ*(5.0/2.0+std::log(2*M_PI*BOLTZ*Temp/(gmx::square(PLANCK))*V/(f*Natom)));
    Shs   = Sig+calc_Shs(f, y);
    rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
    sigHS = std::cbrt(6*y*V/(M_PI*Natom));

    fprintf(fplog, "System = \"%s\"\n", *top.name);
    fprintf(fplog, "Nmol = %d\n", Nmol);
    fprintf(fplog, "Natom = %d\n", Natom);
    fprintf(fplog, "dt = %g ps\n", dt);
    fprintf(fplog, "tmass = %g amu\n", tmass);
    fprintf(fplog, "V = %g nm^3\n", V);
    fprintf(fplog, "rho = %g g/l\n", rho);
    fprintf(fplog, "T = %g K\n", Temp);
    fprintf(fplog, "beta = %g mol/kJ\n", beta);

    fprintf(fplog, "\nDoS parameters\n");
    fprintf(fplog, "Delta = %g\n", Delta);
    fprintf(fplog, "fluidicity = %g\n", f);
    fprintf(fplog, "hard sphere packing fraction = %g\n", y);
    fprintf(fplog, "hard sphere compressibility = %g\n", z);
    fprintf(fplog, "ideal gas entropy = %g\n", Sig);
    fprintf(fplog, "hard sphere entropy = %g\n", Shs);
    fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
    fprintf(fplog, "DoS0 = %g\n", DoS0);
    fprintf(fplog, "Dos2 = %g\n", dos2);
    fprintf(fplog, "DoSTot = %g\n", dostot);

    /* Now compute solid (2) and diffusive (3) components */
    fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
                  bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
                  "\\f{4}S(\\f{12}n\\f{4})", oenv);
    xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
    recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
    for (j = 0; (j < nframes/4); j++)
    {
        dos[DOS_DIFF][j]  = DoS0/(1+gmx::square(DoS0*M_PI*nu[j]/(6*f*Natom)));
        dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
        fprintf(fp, "%10g  %10g  %10g  %10g\n",
                recip_fac*nu[j],
                dos[DOS][j]/recip_fac,
                dos[DOS_SOLID][j]/recip_fac,
                dos[DOS_DIFF][j]/recip_fac);
    }
    xvgrclose(fp);

    /* Finally analyze the results! */
    wCdiff = 0.5;
    wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
    wEdiff = 0.5;
    wAdiff = wEdiff-wSdiff;
    for (j = 0; (j < nframes/4); j++)
    {
        dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
                          dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
        dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +
                          dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
        dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +
                          dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
        dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +
                          dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
    }
    DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
    DiffCoeff = 1000*DiffCoeff/3.0;
    fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
            DiffCoeff);
    fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
            1000*DoS0/(12*tmass*beta));

    cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
                                   nframes/4, &stddev);
    fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
    fprintf(fplog, "\nArrivederci!\n");
    gmx_fio_fclose(fplog);

    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");

    return 0;
}

//.........这里部分代码省略.........
        xvgr_legend(fp, 2, leg, oenv);
        for (i = 0; (i < isize); i++)
        {
            if (!bRes || i+1 == isize ||
                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)
            {
                resind    = top.atoms.atom[index[i]].resind;
                pdb_bfac  = find_pdb_bfac(pdbatoms, &top.atoms.resinfo[resind],
                                          *(top.atoms.atomname[index[i]]));

                fprintf(fp, "%5d  %10.5f  %10.5f\n",
                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, rmsf[i]*bfac,
                        pdb_bfac);
            }
        }
        xvgrclose(fp);
    }
    else
    {
        fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS fluctuation", label, "(nm)", oenv);
        for (i = 0; i < isize; i++)
        {
            if (!bRes || i+1 == isize ||
                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)
            {
                fprintf(fp, "%5d %8.4f\n",
                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, std::sqrt(rmsf[i]));
            }
        }
        xvgrclose(fp);
    }

    for (i = 0; i < isize; i++)
    {
        pdbatoms->pdbinfo[index[i]].bfac = 800*M_PI*M_PI/3.0*rmsf[i];
    }

    if (devfn)
    {
        for (i = 0; i < isize; i++)
        {
            rmsf[i] = (rmsd_x[i][XX]+rmsd_x[i][YY]+rmsd_x[i][ZZ])/count;
        }
        if (bRes)
        {
            average_residues(rmsf, nullptr, 0, isize, index, w_rls, &top.atoms);
        }
        /* Write RMSD output */
        fp = xvgropen(devfn, "RMS Deviation", label, "(nm)", oenv);
        for (i = 0; i < isize; i++)
        {
            if (!bRes || i+1 == isize ||
                top.atoms.atom[index[i]].resind != top.atoms.atom[index[i+1]].resind)
            {
                fprintf(fp, "%5d %8.4f\n",
                        bRes ? top.atoms.resinfo[top.atoms.atom[index[i]].resind].nr : index[i]+1, std::sqrt(rmsf[i]));
            }
        }
        xvgrclose(fp);
    }

    if (opt2bSet("-oq", NFILE, fnm))
    {
        /* Write a .pdb file with B-factors and optionally anisou records */
        for (i = 0; i < isize; i++)
        {
            rvec_inc(pdbx[index[i]], xcm);
        }
        write_sto_conf_indexed(opt2fn("-oq", NFILE, fnm), title, pdbatoms, pdbx,
                               nullptr, ePBC, pdbbox, isize, index);
    }
    if (opt2bSet("-ox", NFILE, fnm))
    {
        rvec *bFactorX;
        snew(bFactorX, top.atoms.nr);
        for (i = 0; i < isize; i++)
        {
            for (d = 0; d < DIM; d++)
            {
                bFactorX[index[i]][d] = xcm[d] + xav[i*DIM + d];
            }
        }
        /* Write a .pdb file with B-factors and optionally anisou records */
        write_sto_conf_indexed(opt2fn("-ox", NFILE, fnm), title, pdbatoms, bFactorX, nullptr,
                               ePBC, pdbbox, isize, index);
        sfree(bFactorX);
    }
    if (bAniso)
    {
        correlate_aniso(opt2fn("-oc", NFILE, fnm), refatoms, pdbatoms, oenv);
        do_view(oenv, opt2fn("-oc", NFILE, fnm), "-nxy");
    }
    do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
    if (devfn)
    {
        do_view(oenv, opt2fn("-od", NFILE, fnm), "-nxy");
    }

    return 0;
}

//.........这里部分代码省略.........
    { efNDX, NULL, NULL,  ffREAD  },
    { efXVG, "-o", "rotacf",  ffWRITE }
  };
#define NFILE asize(fnm)
  int     npargs;
  t_pargs *ppa;
  
  CopyRight(stderr,argv[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,0,NULL);
  
  rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
  
  if (bVec) 
    nvec = isize/2;
  else
    nvec = isize/3;
  
  if (((isize % 3) != 0) && !bVec)
    gmx_fatal(FARGS,"number of index elements not multiple of 3, "
		"these can not be atom triplets\n");
  if (((isize % 2) != 0) && bVec)
    gmx_fatal(FARGS,"number of index elements not multiple of 2, "
		"these can not be atom doublets\n");
  
  top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
  
  snew(c1,nvec);
  for (i=0; (i<nvec); i++)
    c1[i]=NULL;
  n_alloc=0;

  natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
  snew(x_s,natoms);
  
  /* Start the loop over frames */
  t1 = t0 = t;
  teller  = 0;
  do {
    if (teller >= n_alloc) {
      n_alloc+=100;
      for (i=0; (i<nvec); i++)
	srenew(c1[i],DIM*n_alloc);
    }
    t1 = t;
    
    /* Remove periodicity */
    rm_pbc(&(top->idef),ePBC,natoms,box,x,x_s);
    
    /* Compute crossproducts for all vectors, if triplets.
     * else, just get the vectors in case of doublets.
     */
    if (bVec == FALSE) {
      for (i=0; (i<nvec); i++) {
	ai=index[3*i];
	aj=index[3*i+1];
	ak=index[3*i+2];
	rvec_sub(x_s[ai],x_s[aj],xij);
	rvec_sub(x_s[aj],x_s[ak],xjk);
	cprod(xij,xjk,n);
	for(m=0; (m<DIM); m++)
	  c1[i][DIM*teller+m]=n[m];
      }
    }
    else {
      for (i=0; (i<nvec); i++) {
	ai=index[2*i];
	aj=index[2*i+1];
	rvec_sub(x_s[ai],x_s[aj],n);
	for(m=0; (m<DIM); m++)
	  c1[i][DIM*teller+m]=n[m];
      }
    }
    /* Increment loop counter */
    teller++;
  } while (read_next_x(status,&t,natoms,x,box));  
  close_trj(status); 
  fprintf(stderr,"\nDone with trajectory\n");
  
  /* Autocorrelation function */
  if (teller < 2)
    fprintf(stderr,"Not enough frames for correlation function\n");
  else {
    dt=(t1 - t0)/(teller-1);
    
    mode = eacVector;
    
    do_autocorr(ftp2fn(efXVG,NFILE,fnm),"Rotational Correlation Function",
		teller,nvec,c1,dt,mode,bAver);
  }

  do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
    
  thanx(stderr);
    
  return 0;
}

//.........这里部分代码省略.........
                    vol =            (j+1)/(invspz*invspa);
                    break;
                default:
                    vol =             j+1;
                    break;
                }
                if (bMirror)
                    k = j + nradial;
                else
                    k = j;
                grid[i][k] /= nfr*(vol - vol_old);
                if (bMirror)
                    grid[i][nradial-1-j] = grid[i][k];
                vol_old = vol;
                if (grid[i][k] > maxgrid)
                    maxgrid = grid[i][k];
            }
        }
    }
    fprintf(stdout,"\n  The maximum density is %f %s\n",maxgrid,unit);
    if (dmax > 0)
        maxgrid = dmax;

    snew(tickx,n1+1);
    snew(tickz,n2+1);
    if (!bRadial) {
        /* normalize box-axes */
        box1 /= nfr;
        box2 /= nfr;
        for (i=0; i<=n1; i++)
            tickx[i] = i*box1/n1;
        for (i=0; i<=n2; i++)
            tickz[i] = i*box2/n2;
    } else {
        for (i=0; i<=n1; i++)
            tickx[i] = i/invspa - amax;
        if (bMirror) {
            for (i=0; i<=n2; i++)
                tickz[i] = i/invspz - rmax;
        } else {
            for (i=0; i<=n2; i++)
                tickz[i] = i/invspz;
        }
    }

    if (bSums)
    {
        for (i=0; i<n1; ++i)
        {
            fprintf(stdout,"Density sums:\n");
            rowsum=0;
            for (j=0; j<n2; ++j)
                rowsum+=grid[i][j];
            fprintf(stdout,"%g\t",rowsum);
        }
        fprintf(stdout,"\n");
    }

    sprintf(buf,"%s number density",grpname[anagrp]);
    if (!bRadial && (bXmin || bXmax)) {
        if (!bXmax)
            sprintf(buf+strlen(buf),", %c > %g nm",eaver[0][0],xmin);
        else if (!bXmin)
            sprintf(buf+strlen(buf),", %c < %g nm",eaver[0][0],xmax);
        else
            sprintf(buf+strlen(buf),", %c: %g - %g nm",eaver[0][0],xmin,xmax);
    }
    if (ftp2bSet(efDAT,NFILE,fnm))
    {
        fp = ffopen(ftp2fn(efDAT,NFILE,fnm),"w");
        /*optional text form output:  first row is tickz; first col is tickx */
        fprintf(fp,"0\t");
        for(j=0; j<n2; ++j)
            fprintf(fp,"%g\t",tickz[j]);
        fprintf(fp,"\n");

        for (i=0; i<n1; ++i)
        {
            fprintf(fp,"%g\t",tickx[i]);
            for (j=0; j<n2; ++j)
                fprintf(fp,"%g\t",grid[i][j]);
            fprintf(fp,"\n");
        }
        ffclose(fp);
    }
    else
    {
        fp = ffopen(ftp2fn(efXPM,NFILE,fnm),"w");
        write_xpm(fp,MAT_SPATIAL_X | MAT_SPATIAL_Y,buf,unit,
                  bRadial ? "axial (nm)" : label[c1],bRadial ? "r (nm)" : label[c2],
                  n1,n2,tickx,tickz,grid,dmin,maxgrid,rlo,rhi,&nlev);
        ffclose(fp);
    }

    thanx(stderr);

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

    return 0;
}

//.........这里部分代码省略.........
        sum_gyro_tot += sum_gyro;

        if (outp)
        {
            i = -1;
            for (j = 0; j < nat_min/2; j += 2)
            {
                sum_inp[j] /= ninp[j];
                if (i == -1 && sum_inp[j] <= std::exp(-1.0))
                {
                    i = j;
                }
            }
            if (i == -1)
            {
                pers = j;
            }
            else
            {
                /* Do linear interpolation on a log scale */
                pers = i - 2.0
                    + 2.0*(std::log(sum_inp[i-2]) + 1.0)/(std::log(sum_inp[i-2]) - std::log(sum_inp[i]));
            }
            fprintf(outp, "%10.3f %8.4f\n", t*output_env_get_time_factor(oenv), pers);
            sum_pers_tot += pers;
        }

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

    gmx_rmpbc_done(gpbc);

    close_trx(status);

    xvgrclose(out);
    if (outv)
    {
        xvgrclose(outv);
    }
    if (outp)
    {
        xvgrclose(outp);
    }

    sum_eed2_tot /= frame;
    sum_gyro_tot /= frame;
    sum_pers_tot /= frame;
    fprintf(stdout, "\nAverage end to end distance: %.3f (nm)\n",
            std::sqrt(sum_eed2_tot));
    fprintf(stdout, "\nAverage radius of gyration:  %.3f (nm)\n",
            std::sqrt(sum_gyro_tot));
    if (opt2bSet("-p", NFILE, fnm))
    {
        fprintf(stdout, "\nAverage persistence length:  %.2f bonds\n",
                sum_pers_tot);
    }

    /* Handle printing of internal distances. */
    if (outi)
    {
        if (output_env_get_print_xvgr_codes(oenv))
        {
            fprintf(outi, "@    xaxes scale Logarithmic\n");
        }
        ymax = -1;
        ymin = 1e300;
        j    = index[molind[1]-1] - index[molind[0]]; /* Polymer length -1. */
        for (i = 0; i < j; i++)
        {
            intd[i] /= (i + 1) * frame * nmol;
            if (intd[i] > ymax)
            {
                ymax = intd[i];
            }
            if (intd[i] < ymin)
            {
                ymin = intd[i];
            }
        }
        xvgr_world(outi, 1, ymin, j, ymax, oenv);
        for (i = 0; i < j; i++)
        {
            fprintf(outi, "%d  %8.4f\n", i+1, intd[i]);
        }
        xvgrclose(outi);
    }

    do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
    if (opt2bSet("-v", NFILE, fnm))
    {
        do_view(oenv, opt2fn("-v", NFILE, fnm), "-nxy");
    }
    if (opt2bSet("-p", NFILE, fnm))
    {
        do_view(oenv, opt2fn("-p", NFILE, fnm), "-nxy");
    }

    return 0;
}

//.........这里部分代码省略.........
        RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ])),
        norm2(box[YY])==0 ? 0 :
        RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY])));
        printf("new box volume  :%7.2f               (nm^3)\n",det(box));
    }  

    if (check_box(epbcXYZ,box))
        printf("\nWARNING: %s\n",check_box(epbcXYZ,box));

    if (bDist && btype[0][0]=='t')
    {
        if(TRICLINIC(box))
        {
            printf("\nWARNING: Your box is triclinic with non-orthogonal axes. In this case, the\n"
                "distance from the solute to a box surface along the corresponding normal\n"
                "vector might be somewhat smaller than your specified value %f.\n"
                "You can check the actual value with g_mindist -pi\n",dist);
        }
        else
        {
            printf("\nWARNING: No boxtype specified - distance condition applied in each dimension.\n"
                "If the molecule rotates the actual distance will be smaller. You might want\n"
                "to use a cubic box instead, or why not try a dodecahedron today?\n");
        }
    }
    if (bCONECT && (outftp == efPDB) && (inftp == efTPR)) 
        conect = gmx_conect_generate(top);
    else
        conect = NULL;

    if (bIndex) {
        fprintf(stderr,"\nSelect a group for output:\n");
        get_index(&atoms,opt2fn_null("-n",NFILE,fnm),
                  1,&isize,&index,&grpname);
        if (opt2parg_bSet("-label",NPA,pa)) {
            for(i=0; (i<atoms.nr); i++) 
                atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
        }
                
        if (opt2bSet("-bf",NFILE,fnm) || bLegend)
        {
            gmx_fatal(FARGS,"Sorry, cannot do bfactors with an index group.");
        }

        if (outftp == efPDB) 
        {
            out=ffopen(outfile,"w");
            write_pdbfile_indexed(out,title,&atoms,x,ePBC,box,' ',1,isize,index,conect,TRUE);
            ffclose(out);
        }
        else
        {
            write_sto_conf_indexed(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box,isize,index); 
        }
    }
    else {
        if ((outftp == efPDB) || (outftp == efPQR)) {
            out=ffopen(outfile,"w");
            if (bMead) {
                set_pdb_wide_format(TRUE);
                fprintf(out,"REMARK    "
                        "The B-factors in this file hold atomic radii\n");
                fprintf(out,"REMARK    "
                        "The occupancy in this file hold atomic charges\n");
            }
            else if (bGrasp) {
                fprintf(out,"GRASP PDB FILE\nFORMAT NUMBER=1\n");
                fprintf(out,"REMARK    "
                        "The B-factors in this file hold atomic charges\n");
                fprintf(out,"REMARK    "
                        "The occupancy in this file hold atomic radii\n");
            }
            else if (opt2bSet("-bf",NFILE,fnm)) {
                read_bfac(opt2fn("-bf",NFILE,fnm),&n_bfac,&bfac,&bfac_nr);
                set_pdb_conf_bfac(atoms.nr,atoms.nres,&atoms,
                                  n_bfac,bfac,bfac_nr,peratom);
            }
            if (opt2parg_bSet("-label",NPA,pa)) {
                for(i=0; (i<atoms.nr); i++) 
                    atoms.resinfo[atoms.atom[i].resind].chainid=label[0];
            }
            write_pdbfile(out,title,&atoms,x,ePBC,box,' ',-1,conect,TRUE);
            if (bLegend)
                pdb_legend(out,atoms.nr,atoms.nres,&atoms,x);
            if (visbox[0] > 0)
                visualize_box(out,bLegend ? atoms.nr+12 : atoms.nr,
                    bLegend? atoms.nres=12 : atoms.nres,box,visbox);
            ffclose(out);
        }
        else
            write_sto_conf(outfile,title,&atoms,x,bHaveV?v:NULL,ePBC,box); 
    }
    gmx_atomprop_destroy(aps);

    do_view(oenv,outfile,NULL);

    thanx(stderr);

    return 0;
}

//.........这里部分代码省略.........
        { "-mw", FALSE, etBOOL, {&bMW},
          "Use mass weighted fitting" }
    };
    FILE           *out;
    t_trxstatus    *status;
    t_topology      top;
    int             ePBC;
    rvec           *x_ref, *x;
    matrix          box, R;
    real            t;
    int             natoms, i;
    char           *grpname, title[256];
    int             gnx;
    gmx_rmpbc_t     gpbc = NULL;
    atom_id        *index;
    output_env_t    oenv;
    real           *w_rls;
    const char     *leg[]  = { "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" };
#define NLEG asize(leg)
    t_filenm        fnm[] = {
        { efTRX, "-f",   NULL,       ffREAD },
        { efTPS, NULL,   NULL,       ffREAD },
        { efNDX, NULL,   NULL,       ffOPTRD },
        { efXVG, NULL,   "rotmat",   ffWRITE }
    };
#define NFILE asize(fnm)

    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;
    }

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

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

    gmx_rmpbc(gpbc, top.atoms.nr, box, x_ref);

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

    if (reffit[0][0] != 'n')
    {
        get_refx(oenv, ftp2fn(efTRX, NFILE, fnm), reffit[0][2] == 'z' ? 3 : 2, skip,
                 gnx, index, bMW, &top, ePBC, x_ref);
    }

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

    snew(w_rls, natoms);
    for (i = 0; i < gnx; i++)
    {
        if (index[i] >= natoms)
        {
            gmx_fatal(FARGS, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[i]+1, natoms);
        }
        w_rls[index[i]] = (bMW ? top.atoms.atom[index[i]].m : 1.0);
    }

    if (reffit[0][0] == 'n')
    {
        reset_x(gnx, index, natoms, NULL, x_ref, w_rls);
    }

    out = xvgropen(ftp2fn(efXVG, NFILE, fnm),
                   "Fit matrix", "Time (ps)", "", oenv);
    xvgr_legend(out, NLEG, leg, oenv);

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

        reset_x(gnx, index, natoms, NULL, x, w_rls);

        if (bFitXY)
        {
            do_fit_ndim(2, natoms, w_rls, x_ref, x);
        }

        calc_fit_R(DIM, natoms, w_rls, x_ref, x, R);

        fprintf(out,
                "%7g %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n",
                t,
                R[XX][XX], R[XX][YY], R[XX][ZZ],
                R[YY][XX], R[YY][YY], R[YY][ZZ],
                R[ZZ][XX], R[ZZ][YY], R[ZZ][ZZ]);
    }
    while (read_next_x(oenv, status, &t, x, box));

    gmx_rmpbc_done(gpbc);

    close_trj(status);

    gmx_ffclose(out);

    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");

    return 0;
}

int gmx_potential(int argc, char *argv[])
{
    const char        *desc[] = {
        "[THISMODULE] computes the electrostatical potential across the box. The potential is",
        "calculated by first summing the charges per slice and then integrating",
        "twice of this charge distribution. Periodic boundaries are not taken",
        "into account. Reference of potential is taken to be the left side of",
        "the box. It is also possible to calculate the potential in spherical",
        "coordinates as function of r by calculating a charge distribution in",
        "spherical slices and twice integrating them. epsilon_r is taken as 1,",
        "but 2 is more appropriate in many cases."
    };
    gmx_output_env_t  *oenv;
    static int         axis       = 2;       /* normal to memb. default z  */
    static const char *axtitle    = "Z";
    static int         nslices    = 10;      /* nr of slices defined       */
    static int         ngrps      = 1;
    static gmx_bool    bSpherical = FALSE;   /* default is bilayer types   */
    static real        fudge_z    = 0;       /* translate coordinates      */
    static gmx_bool    bCorrect   = 0;
    t_pargs            pa []      = {
        { "-d",   FALSE, etSTR, {&axtitle},
          "Take the normal on the membrane in direction X, Y or Z." },
        { "-sl",  FALSE, etINT, {&nslices},
          "Calculate potential as function of boxlength, dividing the box"
          " in this number of slices." },
        { "-cb",  FALSE, etINT, {&cb},
          "Discard this number of  first slices of box for integration" },
        { "-ce",  FALSE, etINT, {&ce},
          "Discard this number of last slices of box for integration" },
        { "-tz",  FALSE, etREAL, {&fudge_z},
          "Translate all coordinates by this distance in the direction of the box" },
        { "-spherical", FALSE, etBOOL, {&bSpherical},
          "Calculate spherical thingie" },
        { "-ng",       FALSE, etINT, {&ngrps},
          "Number of groups to consider" },
        { "-correct",  FALSE, etBOOL, {&bCorrect},
          "Assume net zero charge of groups to improve accuracy" }
    };
    const char        *bugs[] = {
        "Discarding slices for integration should not be necessary."
    };

    double           **potential,              /* potential per slice        */
    **charge,                                  /* total charge per slice     */
    **field,                                   /* field per slice            */
                       slWidth;                /* width of one slice         */
    char      **grpname;                       /* groupnames                 */
    int        *ngx;                           /* sizes of groups            */
    t_topology *top;                           /* topology        */
    int         ePBC;
    int       **index;                         /* indices for all groups     */
    t_filenm    fnm[] = {                      /* files for g_order       */
        { efTRX, "-f", NULL,  ffREAD },        /* trajectory file             */
        { efNDX, NULL, NULL,  ffREAD },        /* index file          */
        { efTPR, NULL, NULL,  ffREAD },        /* topology file               */
        { efXVG, "-o", "potential", ffWRITE }, /* xvgr output file    */
        { efXVG, "-oc", "charge", ffWRITE },   /* xvgr output file    */
        { efXVG, "-of", "field", ffWRITE },    /* xvgr output file    */
    };

#define NFILE asize(fnm)

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

    /* Calculate axis */
    axis = toupper(axtitle[0]) - 'X';

    top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */

    snew(grpname, ngrps);
    snew(index, ngrps);
    snew(ngx, ngrps);

    rd_index(ftp2fn(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);


    calc_potential(ftp2fn(efTRX, NFILE, fnm), index, ngx,
                   &potential, &charge, &field,
                   &nslices, top, ePBC, axis, ngrps, &slWidth, fudge_z,
                   bSpherical, bCorrect, oenv);

    plot_potential(potential, charge, field, opt2fn("-o", NFILE, fnm),
                   opt2fn("-oc", NFILE, fnm), opt2fn("-of", NFILE, fnm),
                   nslices, ngrps, (const char**)grpname, slWidth, oenv);

    do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);  /* view xvgr file */
    do_view(oenv, opt2fn("-oc", NFILE, fnm), NULL); /* view xvgr file */
    do_view(oenv, opt2fn("-of", NFILE, fnm), NULL); /* view xvgr file */

    return 0;
}

int gmx_h2order(int argc,char *argv[])
{
  const char *desc[] = {
    "Compute the orientation of water molecules with respect to the normal",
    "of the box. The program determines the average cosine of the angle",
    "between de dipole moment of water and an axis of the box. The box is",
    "divided in slices and the average orientation per slice is printed.",
    "Each water molecule is assigned to a slice, per time frame, based on the",
    "position of the oxygen. When -nm  is used the angle between the water",
    "dipole and the axis from the center of mass to the oxygen is calculated",
    "instead of the angle between the dipole and a box axis."
  };
  static int  axis = 2;                       /* normal to memb. default z  */
  static const char *axtitle="Z"; 
  static int  nslices = 0;                    /* nr of slices defined       */
  t_pargs pa[] = {
    { "-d",   FALSE, etSTR, {&axtitle}, 
      "Take the normal on the membrane in direction X, Y or Z." },
    { "-sl",  FALSE, etINT, {&nslices},
      "Calculate order parameter as function of boxlength, dividing the box"
      " in #nr slices."}
  };
  const char *bugs[] = {
    "The program assigns whole water molecules to a slice, based on the first"
    "atom of three in the index file group. It assumes an order O,H,H."
    "Name is not important, but the order is. If this demand is not met,"
    "assigning molecules to slices is different."
  };

  output_env_t oenv;
  real      *slOrder,                       /* av. cosine, per slice      */
            slWidth = 0.0;                  /* width of a slice           */
  rvec      *slDipole;                      
  char      *grpname,            	    /* groupnames                 */
            *micname;
  int       ngx,                            /* nr. of atomsin sol group   */
            nmic;                           /* nr. of atoms in micelle    */
  t_topology *top;                	    /* topology 		  */ 
  int       ePBC;
  atom_id    *index,             	    /* indices for solvent group  */
             *micelle;
  gmx_bool       bMicel =  FALSE;               /* think we're a micel        */
  t_filenm  fnm[] = {             	    /* files for g_order 	  */
    { efTRX, "-f", NULL,  ffREAD },    	    /* trajectory file 	          */
    { efNDX, NULL, NULL,  ffREAD },    	    /* index file 		  */
    { efNDX, "-nm&qu