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

        /***** DIAGONAL ELEMENTS ************/
        for (m = 0; (m < DIM); m++)
        {
            /* Copy the vector data in a linear array */
            for (j = 0; (j < nf2); j++)
            {
                ctmp[j]  = gmx::square(c1[DIM*j+m]);
            }
            if (debug)
            {
                sprintf(buf, "c1diag%d.xvg", m);
                dump_tmp(buf, nf2, ctmp);
            }

            low_do_four_core(nfour, nf2, ctmp, cfour, enNorm);

            if (debug)
            {
                sprintf(buf, "c1dfout%d.xvg", m);
                dump_tmp(buf, nf2, cfour);
            }
            fac = 1.5;
            for (j = 0; (j < nf2); j++)
            {
                csum[j] += fac*(cfour[j]);
            }
        }
        /******* OFF-DIAGONAL ELEMENTS **********/
        for (m = 0; (m < DIM); m++)
        {
            /* Copy the vector data in a linear array */
            m1 = (m+1) % DIM;
            for (j = 0; (j < nf2); j++)
            {
                ctmp[j] = c1[DIM*j+m]*c1[DIM*j+m1];
            }

            if (debug)
            {
                sprintf(buf, "c1off%d.xvg", m);
                dump_tmp(buf, nf2, ctmp);
            }
            low_do_four_core(nfour, nf2, ctmp, cfour, enNorm);
            if (debug)
            {
                sprintf(buf, "c1ofout%d.xvg", m);
                dump_tmp(buf, nf2, cfour);
            }
            fac = 3.0;
            for (j = 0; (j < nf2); j++)
            {
                csum[j] += fac*cfour[j];
            }
        }
    }
    else if (MODE(eacP1) || MODE(eacVector))
    {
        /***************************************************
         * V E C T O R & P1
         ***************************************************/
        if (MODE(eacP1))
        {
            /* First normalize the vectors */
            norm_and_scale_vectors(nframes, c1, 1.0);
        }

        /* For vector thingies we have to do three FFT based correls
         * First for XX, then for YY, then for ZZ
         * After that we sum them and normalise
         */
        for (j = 0; (j < nf2); j++)
        {
            csum[j] = 0.0;
        }
        for (m = 0; (m < DIM); m++)
        {
            /* Copy the vector data in a linear array */
            for (j = 0; (j < nf2); j++)
            {
                ctmp[j] = c1[DIM*j+m];
            }
            low_do_four_core(nfour, nf2, ctmp, cfour, enNorm);
            for (j = 0; (j < nf2); j++)
            {
                csum[j] += cfour[j];
            }
        }
    }
    else
    {
        gmx_fatal(FARGS, "\nUnknown mode in do_autocorr (%d)", mode);
    }

    sfree(cfour);
    for (j = 0; (j < nf2); j++)
    {
        c1[j] = csum[j]/(real)(nframes-j);
    }
}

static void free_params(game_params *params)
{
    sfree(params);
}

static void free_ui(game_ui *ui)
{
    if (ui->sel)
        sfree(ui->sel);
    sfree(ui);
}

static void make_shake_sblock_pd(struct gmx_constr *constr,
				 t_idef *idef,t_mdatoms *md)
{
  int  i,j,m,ncons;
  int  bstart,bnr;
  t_blocka    sblocks;
  t_sortblock *sb;
  t_iatom     *iatom;
  atom_id     *inv_sblock;

  ncons = idef->il[F_CONSTR].nr/3;

  init_blocka(&sblocks);
  gen_sblocks(NULL,md->start,md->start+md->homenr,idef,&sblocks,FALSE);
  
  /*
    bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
    nblocks=blocks->multinr[idef->nodeid] - bstart;
  */
  bstart  = 0;
  constr->nblocks = sblocks.nr;
  if (debug) 
    fprintf(debug,"ncons: %d, bstart: %d, nblocks: %d\n",
	    ncons,bstart,constr->nblocks);
  
  /* Calculate block number for each atom */
  inv_sblock = make_invblocka(&sblocks,md->nr);
  
  done_blocka(&sblocks);
  
  /* Store the block number in temp array and
   * sort the constraints in order of the sblock number 
   * and the atom numbers, really sorting a segment of the array!
   */
#ifdef DEBUGIDEF 
  pr_idef(fplog,0,"Before Sort",idef);
#endif
  iatom=idef->il[F_CONSTR].iatoms;
  snew(sb,ncons);
  for(i=0; (i<ncons); i++,iatom+=3) {
    for(m=0; (m<3); m++)
      sb[i].iatom[m] = iatom[m];
    sb[i].blocknr = inv_sblock[iatom[1]];
  }
  
  /* Now sort the blocks */
  if (debug) {
    pr_sortblock(debug,"Before sorting",ncons,sb);
    fprintf(debug,"Going to sort constraints\n");
  }
  
  qsort(sb,ncons,(size_t)sizeof(*sb),pcomp);
  
  if (debug) {
    fprintf(debug,"I used %d calls to pcomp\n",pcount);
    pr_sortblock(debug,"After sorting",ncons,sb);
  }
  
  iatom=idef->il[F_CONSTR].iatoms;
  for(i=0; (i<ncons); i++,iatom+=3) 
    for(m=0; (m<3); m++)
      iatom[m]=sb[i].iatom[m];
#ifdef DEBUGIDEF
  pr_idef(fplog,0,"After Sort",idef);
#endif
  
  j=0;
  snew(constr->sblock,constr->nblocks+1);
  bnr=-2;
  for(i=0; (i<ncons); i++) {
    if (sb[i].blocknr != bnr) {
      bnr=sb[i].blocknr;
      constr->sblock[j++]=3*i;
    }
  }
  /* Last block... */
  constr->sblock[j++] = 3*ncons;
  
  if (j != (constr->nblocks+1)) {
    fprintf(stderr,"bstart: %d\n",bstart);
    fprintf(stderr,"j: %d, nblocks: %d, ncons: %d\n",
	    j,constr->nblocks,ncons);
    for(i=0; (i<ncons); i++)
      fprintf(stderr,"i: %5d  sb[i].blocknr: %5u\n",i,sb[i].blocknr);
    for(j=0; (j<=constr->nblocks); j++)
      fprintf(stderr,"sblock[%3d]=%5d\n",j,(int)constr->sblock[j]);
    gmx_fatal(FARGS,"DEATH HORROR: "
	      "sblocks does not match idef->il[F_CONSTR]");
  }
  sfree(sb);
  sfree(inv_sblock);
}

void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
                const output_env_t oenv)
{
    FILE        *fp;
    t_mapping   *map;
    int          f, r, *count, *total, ss_count, total_count;
    size_t       s;
    const char** leg;

    map = mat->map;
    snew(count, mat->nmap);
    snew(total, mat->nmap);
    snew(leg, mat->nmap+1);
    leg[0] = "Structure";
    for (s = 0; s < (size_t)mat->nmap; s++)
    {
        leg[s+1] = gmx_strdup(map[s].desc);
    }

    fp = xvgropen(outfile, "Secondary Structure",
                  output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
    if (output_env_get_print_xvgr_codes(oenv))
    {
        fprintf(fp, "@ subtitle \"Structure = ");
    }
    for (s = 0; s < strlen(ss_string); s++)
    {
        if (s > 0)
        {
            fprintf(fp, " + ");
        }
        for (f = 0; f < mat->nmap; f++)
        {
            if (ss_string[s] == map[f].code.c1)
            {
                fprintf(fp, "%s", map[f].desc);
            }
        }
    }
    fprintf(fp, "\"\n");
    xvgr_legend(fp, mat->nmap+1, leg, oenv);

    total_count = 0;
    for (s = 0; s < (size_t)mat->nmap; s++)
    {
        total[s] = 0;
    }
    for (f = 0; f < mat->nx; f++)
    {
        ss_count = 0;
        for (s = 0; s < (size_t)mat->nmap; s++)
        {
            count[s] = 0;
        }
        for (r = 0; r < mat->ny; r++)
        {
            count[mat->matrix[f][r]]++;
            total[mat->matrix[f][r]]++;
        }
        for (s = 0; s < (size_t)mat->nmap; s++)
        {
            if (strchr(ss_string, map[s].code.c1))
            {
                ss_count    += count[s];
                total_count += count[s];
            }
        }
        fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
        for (s = 0; s < (size_t)mat->nmap; s++)
        {
            fprintf(fp, " %5d", count[s]);
        }
        fprintf(fp, "\n");
    }
    /* now print column totals */
    fprintf(fp, "%-8s %5d", "# Totals", total_count);
    for (s = 0; s < (size_t)mat->nmap; s++)
    {
        fprintf(fp, " %5d", total[s]);
    }
    fprintf(fp, "\n");

    /* now print percentages */
    fprintf(fp, "%-8s %5.2f", "# SS %", total_count / (real) (mat->nx * mat->ny));
    for (s = 0; s < (size_t)mat->nmap; s++)
    {
        fprintf(fp, " %5.2f", total[s] / (real) (mat->nx * mat->ny));
    }
    fprintf(fp, "\n");

    gmx_ffclose(fp);
    sfree(leg);
    sfree(count);
}

//.........这里部分代码省略.........
            moltypes[j].res0 = moltypes[j-1].res0+moltypes[j-1].nmol;
        }
        fprintf(stderr, "%7s (%4d atoms): %5d residues\n",
                moltypes[j].name, moltypes[j].natoms, moltypes[j].nmol);
    }

    /* if we have only 1 moleculetype, we don't have to sort */
    if (nrmoltypes > 1)
    {
        /* find out which molecules should go where: */
        moltypes[0].i = moltypes[0].i0 = 0;
        for (j = 1; j < nrmoltypes; j++)
        {
            moltypes[j].i      =
                moltypes[j].i0 =
                    moltypes[j-1].i0+moltypes[j-1].natoms*moltypes[j-1].nmol;
        }

        /* now put them there: */
        snew(newatoms, 1);
        init_t_atoms(newatoms, atoms->nr, FALSE);
        newatoms->nres = atoms->nres;
        snew(newatoms->resinfo, atoms->nres);
        snew(newx, atoms->nr);
        if (v)
        {
            snew(newv, atoms->nr);
        }
        snew(newr, atoms->nr);

        resi_n = 0;
        resnr  = 1;
        j      = 0;
        for (moltp = 0; moltp < nrmoltypes; moltp++)
        {
            i = 0;
            while (i < atoms->nr)
            {
                resi_o = atoms->atom[i].resind;
                if (strcmp(*atoms->resinfo[resi_o].name, moltypes[moltp].name) == 0)
                {
                    /* Copy the residue info */
                    newatoms->resinfo[resi_n]    = atoms->resinfo[resi_o];
                    newatoms->resinfo[resi_n].nr = resnr;
                    /* Copy the atom info */
                    do
                    {
                        newatoms->atom[j]        = atoms->atom[i];
                        newatoms->atomname[j]    = atoms->atomname[i];
                        newatoms->atom[j].resind = resi_n;
                        copy_rvec(x[i], newx[j]);
                        if (v != NULL)
                        {
                            copy_rvec(v[i], newv[j]);
                        }
                        newr[j] = r[i];
                        i++;
                        j++;
                    }
                    while (i < atoms->nr && atoms->atom[i].resind == resi_o);
                    /* Increase the new residue counters */
                    resi_n++;
                    resnr++;
                }
                else
                {
                    /* Skip this residue */
                    do
                    {
                        i++;
                    }
                    while (i < atoms->nr && atoms->atom[i].resind == resi_o);
                }
            }
        }

        /* put them back into the original arrays and throw away temporary arrays */
        sfree(atoms->atomname);
        sfree(atoms->resinfo);
        sfree(atoms->atom);
        sfree(atoms);
        *atoms_solvt = newatoms;
        for (i = 0; i < (*atoms_solvt)->nr; i++)
        {
            copy_rvec(newx[i], x[i]);
            if (v)
            {
                copy_rvec(newv[i], v[i]);
            }
            r[i] = newr[i];
        }
        sfree(newx);
        if (v)
        {
            sfree(newv);
        }
        sfree(newr);
    }
    sfree(moltypes);
}

//.........这里部分代码省略.........
        while (qdelta != 0)
        {
            while (qdelta < 0)
            {
                p_num++;
                qdelta += p_q;
            }
            while (qdelta > 0)
            {
                n_num++;
                qdelta += n_q;
            }
        }
    }

    if ((p_num == 0) && (n_num == 0))
    {
        fprintf(stderr, "No ions to add, will just copy input configuration.\n");
    }
    else
    {
        printf("Will try to add %d %s ions and %d %s ions.\n",
               p_num, p_name, n_num, n_name);
        printf("Select a continuous group of solvent molecules\n");
        get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nwa, &index, &grpname);
        for (i = 1; i < nwa; i++)
        {
            if (index[i] != index[i-1]+1)
            {
                gmx_fatal(FARGS, "The solvent group %s is not continuous: "
                          "index[%d]=%d, index[%d]=%d",
                          grpname, i, index[i-1]+1, i+1, index[i]+1);
            }
        }
        nsa = 1;
        while ((nsa < nwa) &&
               (atoms.atom[index[nsa]].resind ==
                atoms.atom[index[nsa-1]].resind))
        {
            nsa++;
        }
        if (nwa % nsa)
        {
            gmx_fatal(FARGS, "Your solvent group size (%d) is not a multiple of %d",
                      nwa, nsa);
        }
        nw = nwa/nsa;
        fprintf(stderr, "Number of (%d-atomic) solvent molecules: %d\n", nsa, nw);
        if (p_num+n_num > nw)
        {
            gmx_fatal(FARGS, "Not enough solvent for adding ions");
        }

        if (opt2bSet("-p", NFILE, fnm))
        {
            update_topol(opt2fn("-p", NFILE, fnm), p_num, n_num, p_name, n_name, grpname);
        }

        snew(bSet, nw);
        snew(repl, nw);

        snew(v, atoms.nr);
        snew(atoms.pdbinfo, atoms.nr);

        set_pbc(&pbc, ePBC, box);

        if (seed == 0)
        {
            rng = gmx_rng_init(gmx_rng_make_seed());
        }
        else
        {
            rng = gmx_rng_init(seed);
        }
        /* Now loop over the ions that have to be placed */
        while (p_num-- > 0)
        {
            insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
                       1, p_q, p_name, &atoms, rmin, rng);
        }
        while (n_num-- > 0)
        {
            insert_ion(nsa, &nw, bSet, repl, index, x, &pbc,
                       -1, n_q, n_name, &atoms, rmin, rng);
        }
        gmx_rng_destroy(rng);
        fprintf(stderr, "\n");

        if (nw)
        {
            sort_ions(nsa, nw, repl, index, &atoms, x, p_name, n_name);
        }

        sfree(atoms.pdbinfo);
        atoms.pdbinfo = NULL;
    }
    write_sto_conf(ftp2fn(efSTO, NFILE, fnm), *top.name, &atoms, x, NULL, ePBC, box);

    return 0;
}

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

    tmp = realloc(fds, sizeof(*fds) * (fds_num + 1));
    if (tmp == NULL) {
      ERROR("statsd plugin: realloc failed.");
      close(fd);
      continue;
    }
    fds = tmp;
    tmp = fds + fds_num;
    fds_num++;

    memset(tmp, 0, sizeof(*tmp));
    tmp->fd = fd;
    tmp->events = POLLIN | POLLPRI;
    INFO("statsd plugin: Listening on [%s]:%s.", str_node, str_service);
  }

  freeaddrinfo(ai_list);

  if (fds_num == 0) {
    ERROR("statsd plugin: Unable to create listening socket for [%s]:%s.",
          (node != NULL) ? node : "::", service);
    return ENOENT;
  }

  *ret_fds = fds;
  *ret_fds_num = fds_num;
  return 0;
} /* }}} int statsd_network_init */

static void *statsd_network_thread(void *args) /* {{{ */
{
  struct pollfd *fds = NULL;
  size_t fds_num = 0;
  int status;

  status = statsd_network_init(&fds, &fds_num);
  if (status != 0) {
    ERROR("statsd plugin: Unable to open listening sockets.");
    pthread_exit((void *)0);
  }

  while (!network_thread_shutdown) {
    status = poll(fds, (nfds_t)fds_num, /* timeout = */ -1);
    if (status < 0) {

      if ((errno == EINTR) || (errno == EAGAIN))
        continue;

      ERROR("statsd plugin: poll(2) failed: %s", STRERRNO);
      break;
    }

    for (size_t i = 0; i < fds_num; i++) {
      if ((fds[i].revents & (POLLIN | POLLPRI)) == 0)
        continue;

      statsd_network_read(fds[i].fd);
      fds[i].revents = 0;
    }
  } /* while (!network_thread_shutdown) */

  /* Clean up */
  for (size_t i = 0; i < fds_num; i++)
    close(fds[i].fd);
  sfree(fds);

  return (void *)0;
} /* }}} void *statsd_network_thread */

static int statsd_config_timer_percentile(oconfig_item_t *ci) /* {{{ */
{
  double percent = NAN;
  double *tmp;
  int status;

  status = cf_util_get_double(ci, &percent);
  if (status != 0)
    return status;

  if ((percent <= 0.0) || (percent >= 100)) {
    ERROR("statsd plugin: The value for \"%s\" must be between 0 and 100, "
          "exclusively.",
          ci->key);
    return ERANGE;
  }

  tmp =
      realloc(conf_timer_percentile,
              sizeof(*conf_timer_percentile) * (conf_timer_percentile_num + 1));
  if (tmp == NULL) {
    ERROR("statsd plugin: realloc failed.");
    return ENOMEM;
  }
  conf_timer_percentile = tmp;
  conf_timer_percentile[conf_timer_percentile_num] = percent;
  conf_timer_percentile_num++;

  return 0;
} /* }}} int statsd_config_timer_percentile */

//.........这里部分代码省略.........
    bOXT = opt2bSet("-oxt", NFILE, fnm);
    bOV  = opt2bSet("-ov", NFILE, fnm);
    bOF  = opt2bSet("-of", NFILE, fnm);
    bOB  = opt2bSet("-ob", NFILE, fnm);
    bOT  = opt2bSet("-ot", NFILE, fnm);
    bEKT = opt2bSet("-ekt", NFILE, fnm);
    bEKR = opt2bSet("-ekr", NFILE, fnm);
    bCV  = opt2bSet("-cv", NFILE, fnm) || opt2bSet("-av", NFILE, fnm);
    bCF  = opt2bSet("-cf", NFILE, fnm) || opt2bSet("-af", NFILE, fnm);
    bVD  = opt2bSet("-vd", NFILE, fnm) || opt2parg_bSet("-bin", asize(pa), pa);
    if (bMol || bOT || bEKT || bEKR)
    {
        bCom = TRUE;
    }

    bDim[XX]  = bX;
    bDim[YY]  = bY;
    bDim[ZZ]  = bZ;
    bDim[DIM] = bNorm;

    if (bFP)
    {
        sffmt = "\t" gmx_real_fullprecision_pfmt;
    }
    else
    {
        sffmt = "\t%g";
    }
    sprintf(sffmt6, "%s%s%s%s%s%s", sffmt, sffmt, sffmt, sffmt, sffmt, sffmt);

    bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC,
                         &xtop, NULL, topbox,
                         bCom && (bOX || bOXT || bOV || bOT || bEKT || bEKR));
    sfree(xtop);
    if ((bMol || bCV || bCF) && !bTop)
    {
        gmx_fatal(FARGS, "Need a run input file for option -mol, -cv or -cf");
    }

    if (bMol)
    {
        indexfn = ftp2fn(efNDX, NFILE, fnm);
    }
    else
    {
        indexfn = ftp2fn_null(efNDX, NFILE, fnm);
    }

    if (!(bCom && !bMol))
    {
        ngroups = 1;
    }
    snew(grpname, ngroups);
    snew(isize0, ngroups);
    snew(index0, ngroups);
    get_index(&(top.atoms), indexfn, ngroups, isize0, index0, grpname);

    if (bMol)
    {
        mols    = &(top.mols);
        atndx   = mols->index;
        ngroups = isize0[0];
        snew(isize, ngroups);
        snew(index, ngroups);
        for (i = 0; i < ngroups; i++)
        {

/*
 * Helper function to parse a comma-separated list of strings into
 * a preference list array of values. Any missing values are added
 * to the end and duplicates are weeded.
 * XXX: assumes vals in 'mapping' are small +ve integers
 */
static void gprefs(void *sesskey, char *name, char *def,
		   const struct keyvalwhere *mapping, int nvals,
		   Conf *conf, int primary)
{
    char *commalist;
    char *p, *q;
    int i, j, n, v, pos;
    unsigned long seen = 0;	       /* bitmap for weeding dups etc */

    /*
     * Fetch the string which we'll parse as a comma-separated list.
     */
    commalist = gpps_raw(sesskey, name, def);

    /*
     * Go through that list and convert it into values.
     */
    n = 0;
    p = commalist;
    while (1) {
        while (*p && *p == ',') p++;
        if (!*p)
            break;                     /* no more words */

        q = p;
        while (*p && *p != ',') p++;
        if (*p) *p++ = '\0';

        v = key2val(mapping, nvals, q);
        if (v != -1 && !(seen & (1 << v))) {
	    seen |= (1 << v);
            conf_set_int_int(conf, primary, n, v);
            n++;
	}
    }

    sfree(commalist);

    /*
     * Now go through 'mapping' and add values that weren't mentioned
     * in the list we fetched. We may have to loop over it multiple
     * times so that we add values before other values whose default
     * positions depend on them.
     */
    while (n < nvals) {
        for (i = 0; i < nvals; i++) {
	    assert(mapping[i].v < 32);

	    if (!(seen & (1 << mapping[i].v))) {
                /*
                 * This element needs adding. But can we add it yet?
                 */
                if (mapping[i].vrel != -1 && !(seen & (1 << mapping[i].vrel)))
                    continue;          /* nope */

                /*
                 * OK, we can work out where to add this element, so
                 * do so.
                 */
                if (mapping[i].vrel == -1) {
                    pos = (mapping[i].where < 0 ? n : 0);
                } else {
                    for (j = 0; j < n; j++)
                        if (conf_get_int_int(conf, primary, j) ==
                            mapping[i].vrel)
                            break;
                    assert(j < n);     /* implied by (seen & (1<<vrel)) */
                    pos = (mapping[i].where < 0 ? j : j+1);
                }

                /*
                 * And add it.
                 */
                for (j = n-1; j >= pos; j--)
                    conf_set_int_int(conf, primary, j+1,
                                     conf_get_int_int(conf, primary, j));
                conf_set_int_int(conf, primary, pos, mapping[i].v);
                n++;
            }
        }
    }
}

void load_open_settings(void *sesskey, Conf *conf)
{
    int i;
    char *prot;

    conf_set_int(conf, CONF_ssh_subsys, 0);   /* FIXME: load this properly */
    conf_set_str(conf, CONF_remote_cmd, "");
    conf_set_str(conf, CONF_remote_cmd2, "");
    conf_set_str(conf, CONF_ssh_nc_host, "");

    gpps(sesskey, "HostName", "", conf, CONF_host);
    gppfile(sesskey, "LogFileName", conf, CONF_logfilename);
    gppi(sesskey, "LogType", 0, conf, CONF_logtype);
    gppi(sesskey, "LogFileClash", LGXF_ASK, conf, CONF_logxfovr);
    gppi(sesskey, "LogFlush", 1, conf, CONF_logflush);
    gppi(sesskey, "SSHLogOmitPasswords", 1, conf, CONF_logomitpass);
    gppi(sesskey, "SSHLogOmitData", 0, conf, CONF_logomitdata);

    prot = gpps_raw(sesskey, "Protocol", "default");
    conf_set_int(conf, CONF_protocol, default_protocol);
    conf_set_int(conf, CONF_port, default_port);
    {
	const Backend *b = backend_from_name(prot);
	if (b) {
	    conf_set_int(conf, CONF_protocol, b->protocol);
	    gppi(sesskey, "PortNumber", default_port, conf, CONF_port);
	}
    }
    sfree(prot);

    /* Address family selection */
    gppi(sesskey, "AddressFamily", ADDRTYPE_UNSPEC, conf, CONF_addressfamily);

    /* The CloseOnExit numbers are arranged in a different order from
     * the standard FORCE_ON / FORCE_OFF / AUTO. */
    i = gppi_raw(sesskey, "CloseOnExit", 1); conf_set_int(conf, CONF_close_on_exit, (i+1)%3);
    gppi(sesskey, "WarnOnClose", 1, conf, CONF_warn_on_close);
    {
	/* This is two values for backward compatibility with 0.50/0.51 */
	int pingmin, pingsec;
	pingmin = gppi_raw(sesskey, "PingInterval", 0);
	pingsec = gppi_raw(sesskey, "PingIntervalSecs", 0);
	conf_set_int(conf, CONF_ping_interval, pingmin * 60 + pingsec);
    }
    gppi(sesskey, "TCPNoDelay", 1, conf, CONF_tcp_nodelay);
    gppi(sesskey, "TCPKeepalives", 0, conf, CONF_tcp_keepalives);
    gpps(sesskey, "TerminalType", "xterm", conf, CONF_termtype);
    gpps(sesskey, "TerminalSpeed", "38400,38400", conf, CONF_termspeed);
    if (!gppmap(sesskey, "TerminalModes", conf, CONF_ttymodes)) {
	/* This hardcodes a big set of defaults in any new saved
	 * sessions. Let's hope we don't change our mind. */
	for (i = 0; ttymodes[i]; i++)
	    conf_set_str_str(conf, CONF_ttymodes, ttymodes[i], "A");
    }

    /* proxy settings */
    gpps(sesskey, "ProxyExcludeList", "", conf, CONF_proxy_exclude_list);
    i = gppi_raw(sesskey, "ProxyDNS", 1); conf_set_int(conf, CONF_proxy_dns, (i+1)%3);
    gppi(sesskey, "ProxyLocalhost", 0, conf, CONF_even_proxy_localhost);
    gppi(sesskey, "ProxyMethod", -1, conf, CONF_proxy_type);
    if (conf_get_int(conf, CONF_proxy_type) == -1) {
        int i;
        i = gppi_raw(sesskey, "ProxyType", 0);
        if (i == 0)
            conf_set_int(conf, CONF_proxy_type, PROXY_NONE);
        else if (i == 1)
            conf_set_int(conf, CONF_proxy_type, PROXY_HTTP);
        else if (i == 3)
            conf_set_int(conf, CONF_proxy_type, PROXY_TELNET);
        else if (i == 4)
            conf_set_int(conf, CONF_proxy_type, PROXY_CMD);
        else {
            i = gppi_raw(sesskey, "ProxySOCKSVersion", 5);
            if (i == 5)
                conf_set_int(conf, CONF_proxy_type, PROXY_SOCKS5);
            else
                conf_set_int(conf, CONF_proxy_type, PROXY_SOCKS4);
        }
    }
    gpps(sesskey, "ProxyHost", "proxy", conf, CONF_proxy_host);
    gppi(sesskey, "ProxyPort", 80, conf, CONF_proxy_port);
    gpps(sesskey, "ProxyUsername", "", conf, CONF_proxy_username);
    gpps(sesskey, "ProxyPassword", "", conf, CONF_proxy_password);
    gpps(sesskey, "ProxyTelnetCommand", "connect %host %port\\n",
	 conf, CONF_proxy_telnet_command);
    gppmap(sesskey, "Environment", conf, CONF_environmt);
    gpps(sesskey, "UserName", "", conf, CONF_username);
    gpps(sesskey, "Password", "", conf, CONF_password);
    gppi(sesskey, "UserNameFromEnvironment", 0, conf, CONF_username_from_env);
    gpps(sesskey, "LocalUserName", "", conf, CONF_localusername);
    gppi(sesskey, "NoPTY", 0, conf, CONF_nopty);
    gppi(sesskey, "Compression", 0, conf, CONF_compression);
    gppi(sesskey, "TryAgent", 1, conf, CONF_tryagent);
    gppi(sesskey, "AgentFwd", 0, conf, CONF_agentfwd);
    gppi(sesskey, "ChangeUsername", 0, conf, CONF_change_username);
    gppi(sesskey, "GssapiFwd", 0, conf, CONF_gssapifwd);
    gprefs(sesskey, "Cipher", "\0",
	   ciphernames, CIPHER_MAX, conf, CONF_ssh_cipherlist);
    {
	/* Backward-compatibility: we used to have an option to
//.........这里部分代码省略.........

/*
 * Write a set of name/value pairs in the above format, or just the
 * names if include_values is FALSE.
 */
static void wmap(void *handle, char const *outkey, Conf *conf, int primary,
                 int include_values)
{
    char *buf, *p, *q, *key, *realkey, *val;
    int len;

    len = 1;			       /* allow for NUL */

    for (val = conf_get_str_strs(conf, primary, NULL, &key);
	 val != NULL;
	 val = conf_get_str_strs(conf, primary, key, &key))
	len += 2 + 2 * (strlen(key) + strlen(val));   /* allow for escaping */

    buf = snewn(len, char);
    p = buf;

    for (val = conf_get_str_strs(conf, primary, NULL, &key);
	 val != NULL;
	 val = conf_get_str_strs(conf, primary, key, &key)) {

        if (primary == CONF_portfwd && !strcmp(val, "D")) {
            /*
             * Backwards-compatibility hack, as above: translate from
             * the sensible internal representation of dynamic
             * forwardings (key "L<port>", value "D") to the
             * conceptually incoherent legacy storage format (key
             * "D<port>", value empty).
             */
            char *L;

            realkey = key;             /* restore it at end of loop */
            val = "";
            key = dupstr(key);
            L = strchr(key, 'L');
            if (L) *L = 'D';
        } else {
            realkey = NULL;
        }

	if (p != buf)
	    *p++ = ',';
	for (q = key; *q; q++) {
	    if (*q == '=' || *q == ',' || *q == '\\')
		*p++ = '\\';
	    *p++ = *q;
	}
        if (include_values) {
            *p++ = '=';
            for (q = val; *q; q++) {
                if (*q == '=' || *q == ',' || *q == '\\')
                    *p++ = '\\';
                *p++ = *q;
            }
        }

        if (realkey) {
            free(key);
            key = realkey;
        }
    }
    *p = '\0';
    write_setting_s(handle, outkey, buf);
    sfree(buf);
}

/*
 * Read a set of name-value pairs in the format we occasionally use:
 *   NAME\tVALUE\0NAME\tVALUE\0\0 in memory
 *   NAME=VALUE,NAME=VALUE, in storage
 * If there's no "=VALUE" (e.g. just NAME,NAME,NAME) then those keys
 * are mapped to the empty string.
 */
static int gppmap(void *handle, char *name, Conf *conf, int primary)
{
    char *buf, *p, *q, *key, *val;

    /*
     * Start by clearing any existing subkeys of this key from conf.
     */
    while ((key = conf_get_str_nthstrkey(conf, primary, 0)) != NULL)
        conf_del_str_str(conf, primary, key);

    /*
     * Now read a serialised list from the settings and unmarshal it
     * into its components.
     */
    buf = gpps_raw(handle, name, NULL);
    if (!buf)
	return FALSE;

    p = buf;
    while (*p) {
	q = buf;
	val = NULL;
	while (*p && *p != ',') {
	    int c = *p++;
	    if (c == '=')
		c = '\0';
	    if (c == '\\')
		c = *p++;
	    *q++ = c;
	    if (!c)
		val = q;
	}
	if (*p == ',')
	    p++;
	if (!val)
	    val = q;
	*q = '\0';

        if (primary == CONF_portfwd && strchr(buf, 'D') != NULL) {
            /*
             * Backwards-compatibility hack: dynamic forwardings are
             * indexed in the data store as a third type letter in the
             * key, 'D' alongside 'L' and 'R' - but really, they
             * should be filed under 'L' with a special _value_,
             * because local and dynamic forwardings both involve
             * _listening_ on a local port, and are hence mutually
             * exclusive on the same port number. So here we translate
             * the legacy storage format into the sensible internal
             * form, by finding the D and turning it into a L.
             */
            char *newkey = dupstr(buf);
            *strchr(newkey, 'D') = 'L';
            conf_set_str_str(conf, primary, newkey, "D");
            sfree(newkey);
        } else {
            conf_set_str_str(conf, primary, buf, val);
        }
    }
    sfree(buf);

    return TRUE;
}

//.........这里部分代码省略.........
        snew(ctmp, nfour);
    }
    else
    {
        nfour = 0; /* To keep the compiler happy */
        snew(csum, nframes);
        snew(ctmp, nframes);
    }

    /* Loop over items (e.g. molecules or dihedrals)
     * In this loop the actual correlation functions are computed, but without
     * normalizing them.
     */
    for (int i = 0; i < nitem; i++)
    {
        if (bVerbose && (((i % 100) == 0) || (i == nitem-1)))
        {
            fprintf(stderr, "\rThingie %d", i+1);
        }

        if (bFour)
        {
            do_four_core(mode, nfour, nframes, nframes, c1[i], csum, ctmp);
        }
        else
        {
            do_ac_core(nframes, nout, ctmp, c1[i], nrestart, mode);
        }
    }
    if (bVerbose)
    {
        fprintf(stderr, "\n");
    }
    sfree(ctmp);
    sfree(csum);

    if (fn)
    {
        snew(fit, nout);
        fp = xvgropen(fn, title, "Time (ps)", "C(t)", oenv);
    }
    else
    {
        fit = NULL;
        fp  = NULL;
    }
    if (bAver)
    {
        if (nitem > 1)
        {
            average_acf(bVerbose, nframes, nitem, c1);
        }

        if (bNormalize)
        {
            normalize_acf(nout, c1[0]);
        }

        if (eFitFn != effnNONE)
        {
            fit_acf(nout, eFitFn, oenv, fn != NULL, tbeginfit, tendfit, dt, c1[0], fit);
            sum = print_and_integrate(fp, nout, dt, c1[0], fit, 1);
        }
        else
        {
            sum = print_and_integrate(fp, nout, dt, c1[0], NULL, 1);

static void add_solv(const char *fn, t_atoms *atoms, rvec **x, rvec **v, real **exclusionDistances,
                     int ePBC, matrix box,
                     gmx_atomprop_t aps,
                     real defaultDistance, real scaleFactor, int *atoms_added,
                     int *residues_added, real rshell, int max_sol,
                     const output_env_t oenv)
{
    int      i, nmol;
    ivec     n_box;
    char     filename[STRLEN];
    char     title_solvt[STRLEN];
    t_atoms *atoms_solvt;
    rvec    *x_solvt, *v_solvt = NULL;
    real    *exclusionDistances_solvt;
    int      ePBC_solvt;
    matrix   box_solvt;
    int      onr, onres;
    char    *lfn;

    lfn = gmxlibfn(fn);
    strncpy(filename, lfn, STRLEN);
    sfree(lfn);
    {
        int natoms;
        get_stx_coordnum(filename, &natoms);
        if (0 == natoms)
        {
            gmx_fatal(FARGS, "No solvent in %s, please check your input\n", filename);
        }
        snew(atoms_solvt, 1);
        init_t_atoms(atoms_solvt, natoms, FALSE);
    }
    snew(x_solvt, atoms_solvt->nr);
    if (v)
    {
        snew(v_solvt, atoms_solvt->nr);
    }
    snew(exclusionDistances_solvt, atoms_solvt->nr);
    snew(atoms_solvt->resinfo, atoms_solvt->nr);
    snew(atoms_solvt->atomname, atoms_solvt->nr);
    snew(atoms_solvt->atom, atoms_solvt->nr);
    atoms_solvt->pdbinfo = NULL;
    fprintf(stderr, "Reading solvent configuration%s\n",
            v_solvt ? " and velocities" : "");
    read_stx_conf(filename, title_solvt, atoms_solvt, x_solvt, v_solvt,
                  &ePBC_solvt, box_solvt);
    fprintf(stderr, "\"%s\"\n", title_solvt);
    fprintf(stderr, "solvent configuration contains %d atoms in %d residues\n",
            atoms_solvt->nr, atoms_solvt->nres);
    fprintf(stderr, "\n");

    /* apply pbc for solvent configuration for whole molecules */
    rm_res_pbc(atoms_solvt, x_solvt, box_solvt);

    /* initialise distance arrays for solvent configuration */
    exclusionDistances_solvt = makeExclusionDistances(atoms_solvt, aps, defaultDistance, scaleFactor);

    /* calculate the box multiplication factors n_box[0...DIM] */
    nmol = 1;
    for (i = 0; (i < DIM); i++)
    {
        n_box[i] = 1;
        while (n_box[i]*box_solvt[i][i] < box[i][i])
        {
            n_box[i]++;
        }
        nmol *= n_box[i];
    }
    fprintf(stderr, "Will generate new solvent configuration of %dx%dx%d boxes\n",
            n_box[XX], n_box[YY], n_box[ZZ]);

    /* realloc atoms_solvt for the new solvent configuration */
    srenew(atoms_solvt->resinfo, atoms_solvt->nres*nmol);
    srenew(atoms_solvt->atomname, atoms_solvt->nr*nmol);
    srenew(atoms_solvt->atom, atoms_solvt->nr*nmol);
    srenew(x_solvt, atoms_solvt->nr*nmol);
    if (v_solvt)
    {
        srenew(v_solvt, atoms_solvt->nr*nmol);
    }
    srenew(exclusionDistances_solvt, atoms_solvt->nr*nmol);

    /* generate a new solvent configuration */
    make_new_conformation(atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, box_solvt, n_box);

#ifdef DEBUG
    print_stat(x_solvt, atoms_solvt->nr, box_solvt);
#endif

#ifdef DEBUG
    print_stat(x_solvt, atoms_solvt->nr, box_solvt);
#endif
    /* Sort the solvent mixture, not the protein... */
    sort_molecule(&atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt);

    /* add the two configurations */
    onr   = atoms->nr;
    onres = atoms->nres;
    add_conf(atoms, x, v, exclusionDistances, TRUE, ePBC, box, FALSE,
             atoms_solvt, x_solvt, v_solvt, exclusionDistances_solvt, TRUE, rshell, max_sol, oenv);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }
    do
    {
        nframes0++;
        /* calculate x: a fitted struture of the selected atoms */
        if (bPBC)
        {
            gmx_rmpbc(gpbc, nat, box, xread);
        }
        if (bFit)
        {
            reset_x(nfit, ifit, nat, NULL, xread, w_rls);
            do_fit(nat, w_rls, xref, xread);
        }
        for (i = 0; i < natoms; i++)
        {
            rvec_inc(xav[i], xread[index[i]]);
        }
    }
    while (read_next_x(oenv, status, &t, xread, box));
    close_trj(status);

    inv_nframes = 1.0/nframes0;
    for (i = 0; i < natoms; i++)
    {
        for (d = 0; d < DIM; d++)
        {
            xav[i][d]         *= inv_nframes;
            xread[index[i]][d] = xav[i][d];
        }
    }
    write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure",
                           atoms, xread, NULL, epbcNONE, zerobox, natoms, index);
    sfree(xread);

    fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", (int)ndim, (int)ndim);
    nframes = 0;
    nat     = read_first_x(oenv, &status, trxfile, &t, &xread, box);
    tstart  = t;
    do
    {
        nframes++;
        tend = t;
        /* calculate x: a (fitted) structure of the selected atoms */
        if (bPBC)
        {
            gmx_rmpbc(gpbc, nat, box, xread);
        }
        if (bFit)
        {
            reset_x(nfit, ifit, nat, NULL, xread, w_rls);
            do_fit(nat, w_rls, xref, xread);
        }
        if (bRef)
        {
            for (i = 0; i < natoms; i++)
            {
                rvec_sub(xread[index[i]], xref[index[i]], x[i]);
            }
        }
        else
        {
            for (i = 0; i < natoms; i++)
            {
                rvec_sub(xread[index[i]], xav[i], x[i]);
            }

//.........这里部分代码省略.........
    t_pargs         pa[]              = {
        { "-box",    FALSE, etRVEC, {new_box},
          "Box size (in nm)" },
        { "-radius",   FALSE, etREAL, {&defaultDistance},
          "Default van der Waals distance"},
        { "-scale", FALSE, etREAL, {&scaleFactor},
          "Scale factor to multiply Van der Waals radii from the database in share/gromacs/top/vdwradii.dat. The default value of 0.57 yields density close to 1000 g/l for proteins in water." },
        { "-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" },
    };

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

    const char *solventFileName = opt2fn("-cs", NFILE, fnm);
    bProt     = opt2bSet("-cp", NFILE, fnm);
    bBox      = opt2parg_bSet("-box", asize(pa), pa);

    /* check input */
    if (!bProt && !bBox)
    {
        gmx_fatal(FARGS, "When no solute (-cp) is specified, "
                  "a box size (-box) must be specified");
    }

    aps = gmx_atomprop_init();

    snew(atoms, 1);
    init_t_atoms(atoms, 0, FALSE);
    if (bProt)
    {
        /* Generate a solute configuration */
        conf_prot = opt2fn("-cp", NFILE, fnm);
        title     = readConformation(conf_prot, atoms, &x,
                                     bReadV ? &v :