自学教程：C++ srenew函数代码示例

//.........这里部分代码省略.........		    NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);  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++)      y[j][i] = yd[j][i];  }  if (nxtail > nx) {    for(i=nx; (i<nxtail); i++) {      y[0][i] = dt*i+y[0][0];      for(j=1; (j<ny); j++)	y[j][i] = 0.0;    }    nx=nxtail;  }    /* We have read a file WITHOUT standard deviations, so we make our own... */  if (ny==2) {    printf("Creating standard deviation numbers .../n");    srenew(y,3);    snew(y[2],nx);    fac=1.0/((real)nx);    for(i=0; (i<nx); i++)       y[2][i] = fac;  }  eFitFn = sffn2effn(s_ffn);  nfitparm = nfp_ffn[eFitFn];  snew(fitparms,4);  fitparms[0]=tau1;  if (nfitparm > 1)    fitparms[1]=A;  if (nfitparm > 2)    fitparms[2]=tau2;        snew(y[3],nx);  snew(y[4],nx);  snew(y[5],nx);  integral = print_and_integrate(NULL,calc_nbegin(nx,y[0],tbegin),				 dt,y[1],NULL,1);  integral += do_lmfit(nx,y[1],y[2],dt,y[0],tbegin,tend,		       oenv,TRUE,eFitFn,fitparms,fix);  for(i=0; i<nx; i++)    y[3][i] = fit_function(eFitFn,fitparms,y[0][i]);  if (epsRF == 0) {    /* This means infinity! */    lambda = 0;    rffac  = 1;  }  else {    lambda = (eps0 - 1.0)/(2*epsRF - 1.0);    rffac  = (2*epsRF+eps0)/(2*epsRF+1);  }  printf("DATA INTEGRAL: %5.1f, tauD(old) = %5.1f ps, "	  "tau_slope = %5.1f, tau_slope,D = %5.1f ps/n",	  integral,integral*rffac,fitparms[0],fitparms[0]*rffac);  printf("tau_D from tau1 = %8.3g , eps(Infty) = %8.3f/n",	  fitparms[0]*(1 + fitparms[1]*lambda),	  1 + ((1 - fitparms[1])*(eps0 - 1))/(1 + fitparms[1]*lambda));  fitintegral=numerical_deriv(nx,y[0],y[1],y[3],y[4],y[5],tendInt,nsmooth);  printf("FIT INTEGRAL (tau_M): %5.1f, tau_D = %5.1f/n",	  fitintegral,fitintegral*rffac);	    /* Now we have the negative gradient of <Phi(0) Phi(t)> */  write_xvg(opt2fn("-d",NFILE,fnm),"Data",nx-1,6,y,legend,oenv);    /* Do FFT and analysis */    do_four(opt2fn("-o",NFILE,fnm),opt2fn("-c",NFILE,fnm),	  nx-1,y[0],y[5],eps0,epsRF,oenv);  do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");  do_view(oenv,opt2fn("-c",NFILE,fnm),NULL);  do_view(oenv,opt2fn("-d",NFILE,fnm),"-nxy");	      thanx(stderr);  return 0;}

void atoms2md(gmx_mtop_t *mtop, t_inputrec *ir,              int nindex, int *index,              int start, int homenr,              t_mdatoms *md){    gmx_mtop_atomlookup_t alook;    int                   i;    t_grpopts            *opts;    real                  c6, c12;    gmx_groups_t         *groups;    gmx_molblock_t       *molblock;    opts = &ir->opts;    groups = &mtop->groups;    molblock = mtop->molblock;    /* Index==NULL indicates particle decomposition,     * unless we have an empty DD node, so also check for homenr and start.     * This should be signaled properly with an extra parameter or nindex==-1.     */    if (index == NULL && (homenr > 0 || start > 0))    {        md->nr = mtop->natoms;    }    else    {        md->nr = nindex;    }    if (md->nr > md->nalloc)    {        md->nalloc = over_alloc_dd(md->nr);        if (md->nMassPerturbed)        {            srenew(md->massA, md->nalloc);            srenew(md->massB, md->nalloc);        }        srenew(md->massT, md->nalloc);        srenew(md->invmass, md->nalloc);        srenew(md->chargeA, md->nalloc);        srenew(md->c6A, md->nalloc);        srenew(md->sigmaA, md->nalloc);        srenew(md->sigma3A, md->nalloc);        if (md->nPerturbed)        {            srenew(md->chargeB, md->nalloc);            srenew(md->c6B, md->nalloc);            srenew(md->sigmaB, md->nalloc);            srenew(md->sigma3B, md->nalloc);        }        srenew(md->typeA, md->nalloc);        if (md->nPerturbed)        {            srenew(md->typeB, md->nalloc);        }        srenew(md->ptype, md->nalloc);        if (opts->ngtc > 1)        {            srenew(md->cTC, md->nalloc);            /* We always copy cTC with domain decomposition */        }        srenew(md->cENER, md->nalloc);        if (opts->ngacc > 1)        {            srenew(md->cACC, md->nalloc);        }        if (opts->nFreeze &&            (opts->ngfrz > 1 ||             opts->nFreeze[0][XX] || opts->nFreeze[0][YY] || opts->nFreeze[0][ZZ]))        {            srenew(md->cFREEZE, md->nalloc);        }        if (md->bVCMgrps)        {            srenew(md->cVCM, md->nalloc);        }        if (md->bOrires)        {            srenew(md->cORF, md->nalloc);        }        if (md->nPerturbed)        {            srenew(md->bPerturbed, md->nalloc);        }        /* Note that these user t_mdatoms array pointers are NULL         * when there is only one group present.         * Therefore, when adding code, the user should use something like:         * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])         */        if (mtop->groups.grpnr[egcUser1] != NULL)        {            srenew(md->cU1, md->nalloc);        }        if (mtop->groups.grpnr[egcUser2] != NULL)        {            srenew(md->cU2, md->nalloc);//.........这里部分代码省略.........

static void init_lincs(FILE *log,t_topology *top,t_inputrec *ir,		       t_mdatoms *md,int start,int homenr,		       int *nrtot,		       rvec **r,int **bla1,int **bla2,int **blnr,int **blbnb,		       real **bllen,real **blc,real **blcc,real **blm,		       real **tmp1,real **tmp2,real **tmp3,		       real **lincslam,real **bllen0,real **ddist){  t_idef      *idef=&(top->idef);  t_iatom     *iatom;  int         i,j,k,n,b1,b,cen;  int         ncons,nZeroLen;  int         type,a1,a2,b2,nr,n1,n2,nc4;  real        len=0,len1,sign;  real        im1,im2;  int         **at_c,*at_cn,*at_cm;    ncons  = idef->il[F_SHAKE].nr/3;  *nrtot = 0;    if (ncons > 0) {    iatom=idef->il[F_SHAKE].iatoms;    /* Make atom-constraint connection list for temporary use */    snew(at_c,homenr);    snew(at_cn,homenr);    snew(at_cm,homenr);    for(i=0; i<ncons; i++) {      a1=iatom[3*i+1]-start;      a2=iatom[3*i+2]-start;      if (at_cn[a1] >= at_cm[a1]) {	at_cm[a1] += 4;	srenew(at_c[a1],at_cm[a1]);      }      at_c[a1][at_cn[a1]] = i;      at_cn[a1]++;      if (at_cn[a2] >= at_cm[a2]) {	at_cm[a2] += 4;	srenew(at_c[a2],at_cm[a2]);      }      at_c[a2][at_cn[a2]] = i;      at_cn[a2]++;    }    sfree(at_cm);        for(i=0; i<ncons; i++) {      a1=iatom[3*i+1]-start;      a2=iatom[3*i+2]-start;      *nrtot += at_cn[a1] + at_cn[a2] - 2;    }          snew(*r,ncons);    snew(*bla1,ncons);    snew(*bla2,ncons);    snew(*blnr,ncons+1);    snew(*bllen,ncons);    snew(*blc,ncons);    snew(*tmp1,ncons);    snew(*tmp2,ncons);    snew(*tmp3,ncons);    snew(*lincslam,ncons);    snew(*bllen0,ncons);    snew(*ddist,ncons);    snew(*blbnb,*nrtot);    snew(*blcc,*nrtot);    snew(*blm,*nrtot);        /* Make constraint-neighbor list */    (*blnr)[0] = 0;    nZeroLen = 0;    for(i=0; (i<ncons); i++) {      j=3*i;      a1=iatom[j+1];      a2=iatom[j+2];      /* (*blnr)[i+1] = (*blnr)[i] + at_cn[a1] + at_cn[a2] - 2; */      type=iatom[j];      len =idef->iparams[type].shake.dA;      len1=idef->iparams[type].shake.dB;      if (( len<GMX_REAL_MIN || len1<GMX_REAL_MIN) && 	  fabs(len1-len)>GMX_REAL_MIN)	fatal_error(0,"It is not allowed to have a constraint length "		    "zero and non-zero in the A and B topology");      if (len < GMX_REAL_MIN)	nZeroLen++;      (*bla1)[i]=a1;      (*bla2)[i]=a2;      (*bllen)[i]=len;      (*bllen0)[i]=len;      (*ddist)[i]=len1-len;      im1=md->invmass[a1];      im2=md->invmass[a2];      (*blc)[i]=invsqrt(im1+im2);      /* Construct the constraint connection matrix blbnb */      (*blnr)[i+1]=(*blnr)[i];      for(k=0; k<at_cn[a1-start]; k++)	if (at_c[a1-start][k] != i)	  (*blbnb)[((*blnr)[i+1])++]=at_c[a1-start][k];      for(k=0; k<at_cn[a2-start]; k++)//.........这里部分代码省略.........

//.........这里部分代码省略.........        }        if (nparam < 6)        {            fprintf(stderr, "Using default: excluding 3 bonded neighbors/n");            header_settings->nrexcl = 3;        }        if (nparam < 7)        {            fprintf(stderr, "Using default: generating 1,4 H--H interactions/n");            header_settings->bGenerateHH14Interactions = TRUE;        }        if (nparam < 8)        {            fprintf(stderr, "Using default: removing proper dihedrals found on the same bond as a proper dihedral/n");            header_settings->bRemoveDihedralIfWithImproper = TRUE;        }    }    else    {        fprintf(stderr,                "Reading .rtp file without '[ bondedtypes ]' directive,/n"                "Will proceed as if the entry was:/n");        print_resall_header(stderr, header_settings);    }    /* We don't know the current size of rrtp, but simply realloc immediately */    nrtp   = *nrtpptr;    rrtp   = *rtp;    maxrtp = nrtp;    while (!feof(in))    {        if (nrtp >= maxrtp)        {            maxrtp += 100;            srenew(rrtp, maxrtp);        }        /* Initialise rtp entry structure */        rrtp[nrtp] = *header_settings;        if (!get_header(line, header))        {            gmx_fatal(FARGS, "in .rtp file at line:/n%s/n", line);        }        rrtp[nrtp].resname  = gmx_strdup(header);        rrtp[nrtp].filebase = gmx_strdup(filebase);        get_a_line(in, line, STRLEN);        bError       = FALSE;        bNextResidue = FALSE;        do        {            if (!get_header(line, header))            {                bError = TRUE;            }            else            {                bt = get_bt(header);                if (bt != NOTSET)                {                    /* header is an bonded directive */                    bError = !read_bondeds(bt, in, line, &rrtp[nrtp]);                }                else if (gmx_strncasecmp("atoms", header, 5) == 0)                {                    /* header is the atoms directive */                    bError = !read_atoms(in, line, &(rrtp[nrtp]), tab, atype);

static void dielectric(FILE *fmj, FILE *fmd, FILE *outf, FILE *fcur, FILE *mcor,                       FILE *fmjdsp, gmx_bool bNoJump, gmx_bool bACF, gmx_bool bINT,                       int ePBC, t_topology top, t_trxframe fr, real temp,                       real bfit, real efit, real bvit, real evit,                       t_trxstatus *status, int isize, int nmols, int nshift,                       int *index0, int indexm[], real mass2[],                       real qmol[], real eps_rf, const gmx_output_env_t *oenv){    int       i, j;    int       valloc, nalloc, nfr, nvfr;    int       vshfr;    real     *xshfr       = NULL;    int      *vfr         = NULL;    real      refr        = 0.0;    real     *cacf        = NULL;    real     *time        = NULL;    real     *djc         = NULL;    real      corint      = 0.0;    real      prefactorav = 0.0;    real      prefactor   = 0.0;    real      volume;    real      volume_av = 0.0;    real      dk_s, dk_d, dk_f;    real      mj    = 0.0;    real      mj2   = 0.0;    real      mjd   = 0.0;    real      mjdav = 0.0;    real      md2   = 0.0;    real      mdav2 = 0.0;    real      sgk;    rvec      mja_tmp;    rvec      mjd_tmp;    rvec      mdvec;    rvec     *mu    = NULL;    rvec     *xp    = NULL;    rvec     *v0    = NULL;    rvec     *mjdsp = NULL;    real     *dsp2  = NULL;    real      t0;    real      rtmp;    rvec      tmp;    rvec     *mtrans = NULL;    /*     * Variables for the least-squares fit for Einstein-Helfand and Green-Kubo     */    int          bi, ei, ie, ii;    real         rest  = 0.0;    real         sigma = 0.0;    real         malt  = 0.0;    real         err   = 0.0;    real        *xfit;    real        *yfit;    gmx_rmpbc_t  gpbc = NULL;    /*     * indices for EH     */    ei = 0;    bi = 0;    /*     * indices for GK     */    ii  = 0;    ie  = 0;    t0  = 0;    sgk = 0.0;    /* This is the main loop over frames */    nfr = 0;    nvfr   = 0;    vshfr  = 0;    nalloc = 0;    valloc = 0;    clear_rvec(mja_tmp);    clear_rvec(mjd_tmp);    clear_rvec(mdvec);    clear_rvec(tmp);    gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);    do    {        refr = (nfr+1);        if (nfr >= nalloc)        {            nalloc += 100;            srenew(time, nalloc);//.........这里部分代码省略.........

char *wrap_lines(char *buf,int line_width, int indent){  char *b2;  int i,i0,i2,j,b2len,lspace=0,l2space=0;  bool bFirst,bFitsOnLine;  /* characters are copied from buf to b2 with possible spaces changed   * into newlines and extra space added for indentation.   * i indexes buf (source buffer) and i2 indexes b2 (destination buffer)   * i0 points to the beginning of the current line (in buf, source)   * lspace and l2space point to the last space on the current line   * bFirst is set to prevent indentation of first line   * bFitsOnLine says if the first space occurred before line_width, if    * that is not the case, we have a word longer than line_width which    * will also not fit on the next line, so we might as well keep it on    * the current line (where it also won't fit, but looks better)   */    b2=NULL;  b2len=strlen(buf)+1;  snew(b2,b2len);  i0=0;  i2=0;  bFirst=TRUE;  do {    l2space = -1;    /* find the last space before end of line */    for(i=i0; ((i-i0 < line_width) || (l2space==-1)) && (buf[i]); i++) {      b2[i2++] = buf[i];      /* remember the position of a space */      if (buf[i] == ' ') {        lspace = i;	l2space = i2-1;      }      /* if we have a newline before the line is full, reset counters */      if (buf[i]=='/n' && buf[i+1]) { 	i0=i+1;	b2len+=indent;	srenew(b2, b2len);	/* add indentation after the newline */	for(j=0; (j<indent); j++)	  b2[i2++]=' ';      }    }    /* check if one word does not fit on the line */    bFitsOnLine = (i-i0 <= line_width);    /* if we're not at the end of the string */    if (buf[i]) {      /* reset line counters to just after the space */      i0 = lspace+1;      i2 = l2space+1;      /* if the words fit on the line, and we're beyond the indentation part */      if ( (bFitsOnLine) && (l2space >= indent) ) {	/* start a new line */	b2[l2space] = '/n';	/* and add indentation */	if (indent) {	  if (bFirst) {	    line_width-=indent;	    bFirst=FALSE;	  }	  b2len+=indent;	  srenew(b2, b2len);	  for(j=0; (j<indent); j++)	    b2[i2++]=' ';	  /* no extra spaces after indent; */	  while(buf[i0]==' ')	    i0++;	}      }    }  } while (buf[i]);  b2[i2] = '/0';    return b2;}

static void strip_dssp(char *dsspfile,int nres,		       bool bPhobres[],real t,		       real *acc,FILE *fTArea,		       t_matrix *mat,int average_area[]){  static bool bFirst=TRUE;  static char *ssbuf;  FILE *tapeout;  static int xsize,frame;  char buf[STRLEN+1];  char SSTP;  int  i,nr,iacc;  real iaccf,iaccb;  t_xpmelmt c;    tapeout=ffopen(dsspfile,"r");    /* Skip header */  do {    fgets2(buf,STRLEN,tapeout);  } while (strstr(buf,"KAPPA") == NULL);  if (bFirst) {    snew(ssbuf,nres+10);  }    iaccb=iaccf=0;  for(nr=0; (fgets2(buf,STRLEN,tapeout) != NULL); nr++) {    SSTP=buf[16]==' ' ? '~' : buf[16];    ssbuf[nr]=SSTP;        buf[39]='/0';    sscanf(&(buf[34]),"%d",&iacc);    acc[nr]=iacc;    average_area[nr]+=iacc;    if (bPhobres[nr])      iaccb+=iacc;    else      iaccf+=iacc;  }  ssbuf[nr]='/0';    if (bFirst) {    sprintf(mat->title,"Secondary structure");    mat->legend[0]=0;    sprintf(mat->label_x,"%s",time_label());    sprintf(mat->label_y,"Residue");    mat->bDiscrete=TRUE;    mat->ny=nr;    snew(mat->axis_y,nr);    for(i=0; i<nr; i++)      mat->axis_y[i]=i+1;    mat->axis_x=NULL;    mat->matrix=NULL;    xsize=0;    frame=0;    bFirst=FALSE;  }  if (frame>=xsize) {    xsize+=10;    srenew(mat->axis_x,xsize);    srenew(mat->matrix,xsize);  }  mat->axis_x[frame]=t;  snew(mat->matrix[frame],nr);  c.c2=0;  for(i=0; i<nr; i++) {    c.c1=ssbuf[i];    mat->matrix[frame][i] = max(0,searchcmap(mat->nmap,mat->map,c));  }  frame++;  mat->nx=frame;    if (fTArea)    fprintf(fTArea,"%10g  %10g  %10g/n",t,0.01*iaccb,0.01*iaccf);  fclose(tapeout);}

/*! /brief Looks up SETTLE constraints for a range of charge-groups */static void atoms_to_settles(gmx_domdec_t *dd,                             const gmx_mtop_t *mtop,                             const int *cginfo,                             const int **at2settle_mt,                             int cg_start, int cg_end,                             t_ilist *ils_local,                             ind_req_t *ireq){    gmx_ga2la_t            *ga2la;    gmx_mtop_atomlookup_t   alook;    int                     settle;    int                     nral, sa;    int                     cg, a, a_gl, a_glsa, a_gls[3], a_locs[3];    int                     mb, molnr, a_mol, offset;    const gmx_molblock_t   *molb;    const t_iatom          *ia1;    gmx_bool                a_home[3];    int                     nlocal;    gmx_bool                bAssign;    ga2la  = dd->ga2la;    alook = gmx_mtop_atomlookup_settle_init(mtop);    nral = NRAL(F_SETTLE);    for (cg = cg_start; cg < cg_end; cg++)    {        if (GET_CGINFO_SETTLE(cginfo[cg]))        {            for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)            {                a_gl = dd->gatindex[a];                gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);                molb = &mtop->molblock[mb];                settle = at2settle_mt[molb->type][a_mol];                if (settle >= 0)                {                    offset = a_gl - a_mol;                    ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms;                    bAssign = FALSE;                    nlocal  = 0;                    for (sa = 0; sa < nral; sa++)                    {                        a_glsa     = offset + ia1[settle*(1+nral)+1+sa];                        a_gls[sa]  = a_glsa;                        a_home[sa] = ga2la_get_home(ga2la, a_glsa, &a_locs[sa]);                        if (a_home[sa])                        {                            if (nlocal == 0 && a_gl == a_glsa)                            {                                bAssign = TRUE;                            }                            nlocal++;                        }                    }                    if (bAssign)                    {                        if (ils_local->nr+1+nral > ils_local->nalloc)                        {                            ils_local->nalloc = over_alloc_dd(ils_local->nr+1+nral);                            srenew(ils_local->iatoms, ils_local->nalloc);                        }                        ils_local->iatoms[ils_local->nr++] = ia1[settle*4];                        for (sa = 0; sa < nral; sa++)                        {                            if (ga2la_get_home(ga2la, a_gls[sa], &a_locs[sa]))                            {                                ils_local->iatoms[ils_local->nr++] = a_locs[sa];                            }                            else                            {                                ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1;                                /* Add this non-home atom to the list */                                if (ireq->n+1 > ireq->nalloc)                                {                                    ireq->nalloc = over_alloc_large(ireq->n+1);                                    srenew(ireq->ind, ireq->nalloc);                                }                                ireq->ind[ireq->n++] = a_gls[sa];                                /* A check on double atom requests is                                 * not required for settle.                                 */                            }                        }                    }                }            }        }    }//.........这里部分代码省略.........

/*! /brief Looks up constraint for the local atoms */static void atoms_to_constraints(gmx_domdec_t *dd,                                 const gmx_mtop_t *mtop,                                 const int *cginfo,                                 const t_blocka *at2con_mt, int nrec,                                 t_ilist *ilc_local,                                 ind_req_t *ireq){    const t_blocka             *at2con;    gmx_ga2la_t                *ga2la;    gmx_mtop_atomlookup_t       alook;    int                         ncon1;    gmx_molblock_t             *molb;    t_iatom                    *ia1, *ia2, *iap;    int                         nhome, cg, a, a_gl, a_mol, a_loc, b_lo, offset, mb, molnr, b_mol, i, con, con_offset;    gmx_domdec_constraints_t   *dc;    gmx_domdec_specat_comm_t   *dcc;    dc  = dd->constraints;    dcc = dd->constraint_comm;    ga2la  = dd->ga2la;    alook = gmx_mtop_atomlookup_init(mtop);    nhome = 0;    for (cg = 0; cg < dd->ncg_home; cg++)    {        if (GET_CGINFO_CONSTR(cginfo[cg]))        {            for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)            {                a_gl = dd->gatindex[a];                gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);                molb = &mtop->molblock[mb];                ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE);                ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms;                ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms;                /* Calculate the global constraint number offset for the molecule.                 * This is only required for the global index to make sure                 * that we use each constraint only once.                 */                con_offset =                    dc->molb_con_offset[mb] + molnr*dc->molb_ncon_mol[mb];                /* The global atom number offset for this molecule */                offset = a_gl - a_mol;                at2con = &at2con_mt[molb->type];                for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++)                {                    con = at2con->a[i];                    iap = constr_iatomptr(ncon1, ia1, ia2, con);                    if (a_mol == iap[1])                    {                        b_mol = iap[2];                    }                    else                    {                        b_mol = iap[1];                    }                    if (ga2la_get_home(ga2la, offset+b_mol, &a_loc))                    {                        /* Add this fully home constraint at the first atom */                        if (a_mol < b_mol)                        {                            if (dc->ncon+1 > dc->con_nalloc)                            {                                dc->con_nalloc = over_alloc_large(dc->ncon+1);                                srenew(dc->con_gl, dc->con_nalloc);                                srenew(dc->con_nlocat, dc->con_nalloc);                            }                            dc->con_gl[dc->ncon]     = con_offset + con;                            dc->con_nlocat[dc->ncon] = 2;                            if (ilc_local->nr + 3 > ilc_local->nalloc)                            {                                ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3);                                srenew(ilc_local->iatoms, ilc_local->nalloc);                            }                            b_lo = a_loc;                            ilc_local->iatoms[ilc_local->nr++] = iap[0];                            ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a    : b_lo);                            ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a   );                            dc->ncon++;                            nhome++;                        }                    }                    else                    {                        /* We need the nrec constraints coupled to this constraint,                         * so we need to walk out of the home cell by nrec+1 atoms,                         * since already atom bg is not locally present.                         * Therefore we call walk_out with nrec recursions to go                         * after this first call.                         */                        walk_out(con, con_offset, b_mol, offset, nrec,                                 ncon1, ia1, ia2, at2con,//.........这里部分代码省略.........

void connelly_plot(const char *fn,int ndots,real dots[],rvec x[],t_atoms *atoms,		   t_symtab *symtab,int ePBC,matrix box,gmx_bool bSave){  static const char *atomnm="DOT";  static const char *resnm ="DOT";  static const char *title ="Connely Dot Surface Generated by g_sas";  int  i,i0,r0,ii0,k;  rvec *xnew;  t_atoms aaa;  if (bSave) {      i0 = atoms->nr;    r0 = atoms->nres;    srenew(atoms->atom,atoms->nr+ndots);    srenew(atoms->atomname,atoms->nr+ndots);    srenew(atoms->resinfo,r0+1);    atoms->atom[i0].resind = r0;    t_atoms_set_resinfo(atoms,i0,symtab,resnm,r0+1,' ',0,' ');    srenew(atoms->pdbinfo,atoms->nr+ndots);    snew(xnew,atoms->nr+ndots);    for(i=0; (i<atoms->nr); i++)      copy_rvec(x[i],xnew[i]);    for(i=k=0; (i<ndots); i++) {      ii0 = i0+i;      atoms->atomname[ii0] = put_symtab(symtab,atomnm);      atoms->pdbinfo[ii0].type = epdbATOM;      atoms->pdbinfo[ii0].atomnr= ii0;      atoms->atom[ii0].resind = r0;      xnew[ii0][XX] = dots[k++];      xnew[ii0][YY] = dots[k++];      xnew[ii0][ZZ] = dots[k++];      atoms->pdbinfo[ii0].bfac  = 0.0;      atoms->pdbinfo[ii0].occup = 0.0;    }    atoms->nr   = i0+ndots;    atoms->nres = r0+1;    write_sto_conf(fn,title,atoms,xnew,NULL,ePBC,box);    atoms->nres = r0;    atoms->nr   = i0;  }  else {    init_t_atoms(&aaa,ndots,TRUE);    aaa.atom[0].resind = 0;    t_atoms_set_resinfo(&aaa,0,symtab,resnm,1,' ',0,' ');    snew(xnew,ndots);    for(i=k=0; (i<ndots); i++) {      ii0 = i;      aaa.atomname[ii0] = put_symtab(symtab,atomnm);      aaa.pdbinfo[ii0].type = epdbATOM;      aaa.pdbinfo[ii0].atomnr= ii0;      aaa.atom[ii0].resind = 0;      xnew[ii0][XX] = dots[k++];      xnew[ii0][YY] = dots[k++];      xnew[ii0][ZZ] = dots[k++];      aaa.pdbinfo[ii0].bfac  = 0.0;      aaa.pdbinfo[ii0].occup = 0.0;    }    aaa.nr = ndots;    write_sto_conf(fn,title,&aaa,xnew,NULL,ePBC,box);    do_conect(fn,ndots,xnew);    free_t_atoms(&aaa,FALSE);  }  sfree(xnew);}

/*! /brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */static void walk_out(int con, int con_offset, int a, int offset, int nrec,                     int ncon1, const t_iatom *ia1, const t_iatom *ia2,                     const t_blocka *at2con,                     const gmx_ga2la_t *ga2la, gmx_bool bHomeConnect,                     gmx_domdec_constraints_t *dc,                     gmx_domdec_specat_comm_t *dcc,                     t_ilist *il_local,                     ind_req_t *ireq){    int            a1_gl, a2_gl, a_loc, i, coni, b;    const t_iatom *iap;    if (dc->gc_req[con_offset+con] == 0)    {        /* Add this non-home constraint to the list */        if (dc->ncon+1 > dc->con_nalloc)        {            dc->con_nalloc = over_alloc_large(dc->ncon+1);            srenew(dc->con_gl, dc->con_nalloc);            srenew(dc->con_nlocat, dc->con_nalloc);        }        dc->con_gl[dc->ncon]       = con_offset + con;        dc->con_nlocat[dc->ncon]   = (bHomeConnect ? 1 : 0);        dc->gc_req[con_offset+con] = 1;        if (il_local->nr + 3 > il_local->nalloc)        {            il_local->nalloc = over_alloc_dd(il_local->nr+3);            srenew(il_local->iatoms, il_local->nalloc);        }        iap = constr_iatomptr(ncon1, ia1, ia2, con);        il_local->iatoms[il_local->nr++] = iap[0];        a1_gl = offset + iap[1];        a2_gl = offset + iap[2];        /* The following indexing code can probably be optizimed */        if (ga2la_get_home(ga2la, a1_gl, &a_loc))        {            il_local->iatoms[il_local->nr++] = a_loc;        }        else        {            /* We set this index later */            il_local->iatoms[il_local->nr++] = -a1_gl - 1;        }        if (ga2la_get_home(ga2la, a2_gl, &a_loc))        {            il_local->iatoms[il_local->nr++] = a_loc;        }        else        {            /* We set this index later */            il_local->iatoms[il_local->nr++] = -a2_gl - 1;        }        dc->ncon++;    }    /* Check to not ask for the same atom more than once */    if (gmx_hash_get_minone(dc->ga2la, offset+a) == -1)    {        assert(dcc);        /* Add this non-home atom to the list */        if (ireq->n+1 > ireq->nalloc)        {            ireq->nalloc = over_alloc_large(ireq->n+1);            srenew(ireq->ind, ireq->nalloc);        }        ireq->ind[ireq->n++] = offset + a;        /* Temporarily mark with -2, we get the index later */        gmx_hash_set(dc->ga2la, offset+a, -2);    }    if (nrec > 0)    {        for (i = at2con->index[a]; i < at2con->index[a+1]; i++)        {            coni = at2con->a[i];            if (coni != con)            {                /* Walk further */                iap = constr_iatomptr(ncon1, ia1, ia2, coni);                if (a == iap[1])                {                    b = iap[2];                }                else                {                    b = iap[1];                }                if (!ga2la_get_home(ga2la, offset+b, &a_loc))                {                    walk_out(coni, con_offset, b, offset, nrec-1,                             ncon1, ia1, ia2, at2con,                             ga2la, FALSE, dc, dcc, il_local, ireq);                }            }        }    }}

//.........这里部分代码省略.........    printf("Root mean square deviation after lsq fit = %g nm/n", rms);    if (bBfac)    {        printf("Atomic MSD's range from %g to %g nm^2/n", minmsd, maxmsd);    }    if (bFit)    {        /* reset coordinates of reference and fitted structure */        for (i = 0; i < atoms1->nr; i++)        {            for (m = 0; m < DIM; m++)            {                x1[i][m] += xcm1[m];            }        }        for (i = 0; i < atoms2->nr; i++)        {            for (m = 0; m < DIM; m++)            {                x2[i][m] += xcm1[m];            }        }    }    outfile = ftp2fn(efSTO, NFILE, fnm);    switch (fn2ftp(outfile))    {        case efPDB:        case efBRK:        case efENT:            if (bBfac || bLabel)            {                srenew(atoms1->pdbinfo, atoms1->nr);                srenew(atoms1->atom, atoms1->nr); /* Why renew atom? */                /* Avoid segfaults when writing the pdb-file */                for (i = 0; i < atoms1->nr; i++)                {                    atoms1->pdbinfo[i].type         = eptAtom;                    atoms1->pdbinfo[i].occup        = 1.00;                    atoms1->pdbinfo[i].bAnisotropic = FALSE;                    if (bBfac)                    {                        atoms1->pdbinfo[i].bfac = 0;                    }                    if (bLabel)                    {                        atoms1->resinfo[atoms1->atom[i].resind].chainid = 'A';                    }                }                for (i = 0; i < isize1; i++)                {                    /* atoms1->pdbinfo[index1[i]].type = eptAtom; *//*  atoms1->pdbinfo[index1[i]].bAnisotropic = FALSE; */                    if (bBfac)                    {                        atoms1->pdbinfo[index1[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];                    }/*  if (bLabel) *//*    atoms1->resinfo[atoms1->atom[index1[i]].resind].chain = 'A'; */                }                srenew(atoms2->pdbinfo, atoms2->nr);                srenew(atoms2->atom, atoms2->nr); /* Why renew atom? */

void set_lincs(t_idef *idef,t_mdatoms *md,               gmx_bool bDynamics,t_commrec *cr,               struct gmx_lincsdata *li){    int      start,natoms,nflexcon;    t_blocka at2con;    t_iatom  *iatom;    int      i,k,ncc_alloc,ni,con,nconnect,concon;    int      type,a1,a2;    real     lenA=0,lenB;    gmx_bool     bLocal;    li->nc = 0;    li->ncc = 0;		    /* This is the local topology, so there are only F_CONSTR constraints */    if (idef->il[F_CONSTR].nr == 0)    {        /* There are no constraints,         * we do not need to fill any data structures.         */        return;    }        if (debug)    {        fprintf(debug,"Building the LINCS connectivity/n");    }        if (DOMAINDECOMP(cr))    {        if (cr->dd->constraints)        {            dd_get_constraint_range(cr->dd,&start,&natoms);        }        else        {            natoms = cr->dd->nat_home;        }        start = 0;    }    else if(PARTDECOMP(cr))	{		pd_get_constraint_range(cr->pd,&start,&natoms);	}	else    {        start  = md->start;        natoms = md->homenr;    }    at2con = make_at2con(start,natoms,idef->il,idef->iparams,bDynamics,                         &nflexcon);	    if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0)    {        li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3);        srenew(li->bllen0,li->nc_alloc);        srenew(li->ddist,li->nc_alloc);        srenew(li->bla,2*li->nc_alloc);        srenew(li->blc,li->nc_alloc);        srenew(li->blc1,li->nc_alloc);        srenew(li->blnr,li->nc_alloc+1);        srenew(li->bllen,li->nc_alloc);        srenew(li->tmpv,li->nc_alloc);        srenew(li->tmp1,li->nc_alloc);        srenew(li->tmp2,li->nc_alloc);        srenew(li->tmp3,li->nc_alloc);        srenew(li->lambda,li->nc_alloc);        if (li->ncg_triangle > 0)        {            /* This is allocating too much, but it is difficult to improve */            srenew(li->triangle,li->nc_alloc);            srenew(li->tri_bits,li->nc_alloc);        }    }        iatom = idef->il[F_CONSTR].iatoms;        ncc_alloc = li->ncc_alloc;    li->blnr[0] = 0;        ni = idef->il[F_CONSTR].nr/3;    con = 0;    nconnect = 0;    li->blnr[con] = nconnect;    for(i=0; i<ni; i++)    {        bLocal = TRUE;        type = iatom[3*i];        a1   = iatom[3*i+1];        a2   = iatom[3*i+2];        lenA = idef->iparams[type].constr.dA;        lenB = idef->iparams[type].constr.dB;        /* Skip the flexible constraints when not doing dynamics */        if (bDynamics || lenA!=0 || lenB!=0)        {            li->bllen0[con]  = lenA;            li->ddist[con]   = lenB - lenA;//.........这里部分代码省略.........

//.........这里部分代码省略.........                ii = pref->params.ind[i];                if (ga2la)                {                    if (!ga2la_get_home(ga2la, pref->params.ind[i], &ii))                    {                        ii = -1;                    }                }                if (ii >= start && ii < end)                {                    double dr2, dr2_rel, inp;                    dvec   dr;                    pbc_dx_aiuc(pbc, x[ii], g_x, dx);                    inp = iprod(dir, dx);                    dr2 = 0;                    for (m = 0; m < DIM; m++)                    {                        /* Determine the radial components */                        dr[m] = dx[m] - inp*dir[m];                        dr2  += dr[m]*dr[m];                    }                    dr2_rel = dr2*inv_cyl_r2;                    if (dr2_rel < 1)                    {                        double mass, weight, dweight_r;                        dvec   mdw;                        /* add to index, to sum of COM, to weight array */                        if (pdyna->nat_loc >= pdyna->nalloc_loc)                        {                            pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);                            srenew(pdyna->ind_loc,    pdyna->nalloc_loc);                            srenew(pdyna->weight_loc, pdyna->nalloc_loc);                            srenew(pdyna->mdw,        pdyna->nalloc_loc);                            srenew(pdyna->dv,         pdyna->nalloc_loc);                        }                        pdyna->ind_loc[pdyna->nat_loc] = ii;                        mass      = md->massT[ii];                        /* The radial weight function is 1-2x^2+x^4,                         * where x=r/cylinder_r. Since this function depends                         * on the radial component, we also get radial forces                         * on both groups.                         */                        weight    = 1 + (-2 + dr2_rel)*dr2_rel;                        dweight_r = (-4 + 4*dr2_rel)*inv_cyl_r2;                        pdyna->weight_loc[pdyna->nat_loc] = weight;                        sum_a    += mass*weight*inp;                        wmass    += mass*weight;                        wwmass   += mass*weight*weight;                        dsvmul(mass*dweight_r, dr, mdw);                        copy_dvec(mdw, pdyna->mdw[pdyna->nat_loc]);                        /* Currently we only have the axial component of the                         * distance (inp) up to an unkown offset. We add this                         * offset after the reduction needs to determine the                         * COM of the cylinder group.                         */                        pdyna->dv[pdyna->nat_loc] = inp;                        for (m = 0; m < DIM; m++)                        {                            radf_fac0[m] += mdw[m];                            radf_fac1[m] += mdw[m]*inp;                        }                        pdyna->nat_loc++;

static void atomcat(t_atoms *dest, t_atoms *src, int copies,                    int maxres_renum, int *maxresnr){    int i,j,l,size;    int srcnr=src->nr;    int destnr=dest->nr;    if (srcnr)    {        size=destnr+copies*srcnr;        srenew(dest->atom,size);        srenew(dest->atomname,size);        srenew(dest->atomtype,size);        srenew(dest->atomtypeB,size);    }    if (src->nres)    {        size=dest->nres+copies*src->nres;        srenew(dest->resinfo,size);    }        /* residue information */    for (l=dest->nres,j=0; (j<copies); j++,l+=src->nres)    {        memcpy((char *) &(dest->resinfo[l]),(char *) &(src->resinfo[0]),               (size_t)(src->nres*sizeof(src->resinfo[0])));    }        for (l=destnr,j=0; (j<copies); j++,l+=srcnr)    {        memcpy((char *) &(dest->atomname[l]),(char *) &(src->atomname[0]),               (size_t)(srcnr*sizeof(src->atomname[0])));        memcpy((char *) &(dest->atomtype[l]),(char *) &(src->atomtype[0]),               (size_t)(srcnr*sizeof(src->atomtype[0])));        memcpy((char *) &(dest->atomtypeB[l]),(char *) &(src->atomtypeB[0]),               (size_t)(srcnr*sizeof(src->atomtypeB[0])));        memcpy((char *) &(dest->atom[l]),(char *) &(src->atom[0]),               (size_t)(srcnr*sizeof(src->atom[0])));    }        /* Increment residue indices */    for (l=destnr,j=0; (j<copies); j++)    {        for (i=0; (i<srcnr); i++,l++)        {            dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;        }    }            if (src->nres <= maxres_renum)    {        /* Single residue molecule, continue counting residues */        for (j=0; (j<copies); j++)        {            for (l=0; l<src->nres; l++)            {                (*maxresnr)++;                dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;            }        }    }        dest->nres += copies*src->nres;    dest->nr   += copies*src->nr;}

//.........这里部分代码省略.........                    else                    {                        ilst = &dc->ils[thread];                    }                    ilst->nr = 0;                    ireqt = &dd->constraint_comm->ireq[thread];                    if (thread > 0)                    {                        ireqt->n = 0;                    }                    atoms_to_settles(dd, mtop, cginfo, at2settle_mt,                                     cg0, cg1,                                     ilst, ireqt);                }            }            GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;        }        /* Combine the generate settles and requested indices */        for (thread = 1; thread < dc->nthread; thread++)        {            t_ilist   *ilst;            ind_req_t *ireqt;            int        ia;            if (thread > t0_set)            {                ilst = &dc->ils[thread];                if (ils_local->nr + ilst->nr > ils_local->nalloc)                {                    ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr);                    srenew(ils_local->iatoms, ils_local->nalloc);                }                for (ia = 0; ia < ilst->nr; ia++)                {                    ils_local->iatoms[ils_local->nr+ia] = ilst->iatoms[ia];                }                ils_local->nr += ilst->nr;            }            ireqt = &dd->constraint_comm->ireq[thread];            if (ireq->n+ireqt->n > ireq->nalloc)            {                ireq->nalloc = over_alloc_large(ireq->n+ireqt->n);                srenew(ireq->ind, ireq->nalloc);            }            for (ia = 0; ia < ireqt->n; ia++)            {                ireq->ind[ireq->n+ia] = ireqt->ind[ia];            }            ireq->n += ireqt->n;        }        if (debug)        {            fprintf(debug, "Settles: total %3d/n", ils_local->nr/4);        }    }    if (dd->constraint_comm)    {        int nral1;        at_end =

//.........这里部分代码省略.........                }            }        }    }    fprintf(stderr, "/n");    nset = n;    snew(eneset, nset+1);    fprintf(stderr, "Will select half-matrix of energies with %d elements/n", n);    /* Start reading energy frames */    snew(fr, 1);    do    {        do        {            bCont = do_enx(in, fr);            if (bCont)            {                timecheck = check_times(fr->t);            }        }        while (bCont && (timecheck < 0));        if (timecheck == 0)        {#define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE            if (bCont)            {                fprintf(stderr, "/rRead frame: %d, Time: %.3f", teller, fr->t);                if ((nenergy % 1000) == 0)                {                    srenew(time, nenergy+1000);                    for (i = 0; (i <= nset); i++)                    {                        srenew(eneset[i], nenergy+1000);                    }                }                time[nenergy] = fr->t;                sum           = 0;                for (i = 0; (i < nset); i++)                {                    eneset[i][nenergy] = fr->ener[set[i]].e;                    sum               += fr->ener[set[i]].e;                }                if (bSum)                {                    eneset[nset][nenergy] = sum;                }                nenergy++;            }            teller++;        }    }    while (bCont && (timecheck == 0));    fprintf(stderr, "/n");    fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "            "over %d frames/n", ngroups, nset, nenergy);    snew(emat, egNR+egSP);    for (j = 0; (j < egNR+egSP); j++)    {        if (egrp_use[m])

//.........这里部分代码省略.........    if (opt2parg_bSet("-rt", NPA, pa))    {        otfile  = opt2fn("-ot", NFILE, fnm);    }    else    {        otfile  = opt2fn_null("-ot", NFILE, fnm);    }    if (!matfile && !otfile && !orfile)    {        fprintf(stderr,                "For output set one (or more) of the output file options/n");        exit(0);    }    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, boxtop,                  FALSE);    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);    nalloc = 0;    time   = NULL;    sbox   = NULL;    sx     = NULL;    clear_mat(avbox);    natom = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);    nfr   = 0;    do    {        if (nfr >= nalloc)        {            nalloc += 100;            srenew(time, nalloc);            srenew(sbox, nalloc);            srenew(sx, nalloc);        }        assert(time != NULL); assert(sbox != NULL);        time[nfr] = t;        copy_mat(box, sbox[nfr]);        /* This assumes that the off-diagonal box elements         * are not affected by jumps across the periodic boundaries.         */        m_add(avbox, box, avbox);        snew(sx[nfr], isize);        for (i = 0; i < isize; i++)        {            copy_rvec(x[index[i]], sx[nfr][i]);        }        nfr++;    }    while (read_next_x(oenv, status, &t, x, box));    /* clean up */    sfree(x);    close_trj(status);    fprintf(stderr, "Read %d frames/n", nfr);    dt = (time[nfr-1] - time[0])/(nfr - 1);    /* Some ugly rounding to get nice nice times in the output */    dt = (int)(10000.0*dt + 0.5)/10000.0;    invbin = 1.0/rbin;

//.........这里部分代码省略.........    gmx_fatal(FARGS,"DSSP executable (%s) does not exist (use setenv DSSP)",		dptr);  sprintf(dssp,"%s %s %s %s > /dev/null %s",	  dptr,bDoAccSurf?"":"-na",pdbfile,tmpfile,bVerbose?"":"2> /dev/null");  if (bVerbose)    fprintf(stderr,"dssp cmd='%s'/n",dssp);    if (fnTArea) {    fTArea=xvgropen(fnTArea,"Solvent Accessible Surface Area",		    xvgr_tlabel(),"Area (nm//S2//N)");    xvgr_legend(fTArea,2,leg);  } else    fTArea=NULL;    mat.map=NULL;  mat.nmap=getcmap(libopen(opt2fn("-map",NFILE,fnm)),		   opt2fn("-map",NFILE,fnm),&(mat.map));    natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);  if (natoms > atoms->nr)     gmx_fatal(FARGS,"/nTrajectory does not match topology!");  if (gnx > natoms)    gmx_fatal(FARGS,"/nTrajectory does not match selected group!");    snew(average_area,atoms->nres+10);  snew(av_area,atoms->nres+10);  snew(norm_av_area,atoms->nres+10);  accr=NULL;  naccr=0;  do {    t = convert_time(t);    if (nframe>=naccr) {      naccr+=10;      srenew(accr,naccr);      for(i=naccr-10; i<naccr; i++)	snew(accr[i],atoms->nres+10);    }    rm_pbc(&(top.idef),ePBC,natoms,box,x,x);    tapein=ffopen(pdbfile,"w");    write_pdbfile_indexed(tapein,NULL,atoms,x,ePBC,box,0,-1,gnx,index);    fclose(tapein);#ifdef GMX_NO_SYSTEM    printf("Warning-- No calls to system(3) supported on this platform.");    printf("Warning-- Skipping execution of 'system(/"%s/")'.", dssp);    exit(1);#else    if(0 != system(dssp))    {	gmx_fatal(FARGS,"Failed to execute command: %s",dssp);    }#endif    strip_dssp(tmpfile,nres,bPhbres,t,	       accr[nframe],fTArea,&mat,average_area);    remove(tmpfile);    remove(pdbfile);    nframe++;  } while(read_next_x(status,&t,natoms,x,box));  fprintf(stderr,"/n");  close_trj(status);  if (fTArea)    ffclose(fTArea);    prune_ss_legend(&mat);  

//.........这里部分代码省略.........        {            string += 8;            if (check_have_atoms(atoms, ostring) &&                parse_int(&string, &sel_nr) &&                (sel_nr >= 0) && (sel_nr < block->nr))            {                split_group(atoms, sel_nr, block, gn, FALSE);            }        }        else if (strncmp(string, "splitat", 7) == 0)        {            string += 7;            if (check_have_atoms(atoms, ostring) &&                parse_int(&string, &sel_nr) &&                (sel_nr >= 0) && (sel_nr < block->nr))            {                split_group(atoms, sel_nr, block, gn, TRUE);            }        }        else if (string[0] == '/0')        {            bPrintOnce = TRUE;        }        else if (string[0] != 'q')        {            nr1 = -1;            nr2 = -1;            if (parse_entry(&string, natoms, atoms, block, gn, &nr, index, gname))            {                do                {                    while (string[0] == ' ')                    {                        string++;                    }                    bAnd = FALSE;                    bOr  = FALSE;                    if (string[0] == '&')                    {                        bAnd = TRUE;                    }                    else if (string[0] == '|')                    {                        bOr = TRUE;                    }                    if (bAnd || bOr)                    {                        string++;                        nr1 = nr;                        for (i = 0; i < nr; i++)                        {                            index1[i] = index[i];                        }                        strcpy(gname1, gname);                        if (parse_entry(&string, natoms, atoms, block, gn, &nr2, index2, gname2))                        {                            if (bOr)                            {                                or_groups(nr1, index1, nr2, index2, &nr, index);                                sprintf(gname, "%s_%s", gname1, gname2);                            }                            else                            {                                and_groups(nr1, index1, nr2, index2, &nr, index);                                sprintf(gname, "%s_&_%s", gname1, gname2);                            }                        }                    }                }                while (bAnd || bOr);            }            while (string[0] == ' ')            {                string++;            }            if (string[0])            {                printf("/nSyntax error: /"%s/"/n", string);            }            else if (nr > 0)            {                copy2block(nr, index, block);                srenew(*gn, block->nr);                newgroup        = block->nr-1;                (*gn)[newgroup] = strdup(gname);            }            else            {                printf("Group is empty/n");            }        }    }    while (string[0] != 'q');    sfree(index);    sfree(index1);    sfree(index2);}

int gmx_saltbr(int argc, char *argv[]){    const char     *desc[] = {        "[THISMODULE] plots the distance between all combination of charged groups",        "as a function of time. The groups are combined in different ways.",        "A minimum distance can be given (i.e. a cut-off), such that groups",        "that are never closer than that distance will not be plotted.[PAR]",        "Output will be in a number of fixed filenames, [TT]min-min.xvg[tt], [TT]plus-min.xvg[tt]",        "and [TT]plus-plus.xvg[tt], or files for every individual ion pair if the [TT]-sep[tt]",        "option is selected. In this case, files are named as [TT]sb-(Resname)(Resnr)-(Atomnr)[tt].",        "There may be [BB]many[bb] such files."    };    static gmx_bool bSep     = FALSE;    static real     truncate = 1000.0;    t_pargs         pa[]     = {        { "-t",   FALSE, etREAL, {&truncate},          "Groups that are never closer than this distance are not plotted" },        { "-sep", FALSE, etBOOL, {&bSep},          "Use separate files for each interaction (may be MANY)" }    };    t_filenm        fnm[] = {        { efTRX, "-f",  NULL, ffREAD },        { efTPR, NULL,  NULL, ffREAD },    };#define NFILE asize(fnm)    FILE              *out[3], *fp;    static const char *title[3] = {        "Distance between positively charged groups",        "Distance between negatively charged groups",        "Distance between oppositely charged groups"    };    static const char *fn[3] = {        "plus-plus.xvg",        "min-min.xvg",        "plus-min.xvg"    };    int                nset[3] = {0, 0, 0};    t_topology        *top;    int                ePBC;    char              *buf;    t_trxstatus       *status;    int                i, j, k, m, nnn, teller, ncg;    real               t, *time, qi, qj;    t_charge          *cg;    real            ***cgdist;    int              **nWithin;    t_pbc              pbc;    rvec              *x;    matrix             box;    gmx_output_env_t  *oenv;    if (!parse_common_args(&argc, argv, PCA_CAN_TIME,                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))    {        return 0;    }    top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC);    cg  = mk_charge(&top->atoms, &(top->cgs), &ncg);    snew(cgdist, ncg);    snew(nWithin, ncg);    for (i = 0; (i < ncg); i++)    {        snew(cgdist[i], ncg);        snew(nWithin[i], ncg);    }    read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);    teller = 0;    time   = NULL;    do    {        srenew(time, teller+1);        time[teller] = t;        set_pbc(&pbc, ePBC, box);        for (i = 0; (i < ncg); i++)        {            for (j = i+1; (j < ncg); j++)            {                srenew(cgdist[i][j], teller+1);                cgdist[i][j][teller] =                    calc_dist(&pbc, x, &(top->cgs), cg[i].cg, cg[j].cg);                if (cgdist[i][j][teller] < truncate)                {                    nWithin[i][j] = 1;                }            }        }        teller++;    }    while (read_next_x(oenv, status, &t, x, box));    fprintf(stderr, "/n");    close_trj(status);//.........这里部分代码省略.........

//.........这里部分代码省略.........    int          nndxin;    const char  *stxfile;    char       **ndxinfiles;    const char  *ndxoutfile;    gmx_bool     bNatoms;    int          i, j;    t_atoms     *atoms;    rvec        *x, *v;    int          ePBC;    matrix       box;    t_blocka    *block, *block2;    char       **gnames, **gnames2;    t_filenm     fnm[] = {        { efSTX, "-f", NULL,     ffOPTRD  },        { efNDX, "-n", NULL,     ffOPTRDMULT },        { efNDX, "-o", NULL,     ffWRITE }    };#define NFILE asize(fnm)    if (!parse_common_args(&argc, argv, 0, NFILE, fnm, NPA, pa, asize(desc), desc,                           0, NULL, &oenv))    {        return 0;    }    stxfile = ftp2fn_null(efSTX, NFILE, fnm);    if (opt2bSet("-n", NFILE, fnm))    {        nndxin = opt2fns(&ndxinfiles, "-n", NFILE, fnm);    }    else    {        nndxin = 0;    }    ndxoutfile = opt2fn("-o", NFILE, fnm);    bNatoms    = opt2parg_bSet("-natoms", NPA, pa);    if (!stxfile && !nndxin)    {        gmx_fatal(FARGS, "No input files (structure or index)");    }    if (stxfile)    {        snew(atoms, 1);        get_stx_coordnum(stxfile, &(atoms->nr));        init_t_atoms(atoms, atoms->nr, TRUE);        snew(x, atoms->nr);        snew(v, atoms->nr);        fprintf(stderr, "/nReading structure file/n");        read_stx_conf(stxfile, title, atoms, x, v, &ePBC, box);        natoms  = atoms->nr;        bNatoms = TRUE;    }    else    {        atoms = NULL;        x     = NULL;    }    /* read input file(s) */    block  = new_blocka();    gnames = NULL;    printf("Going to read %d old index file(s)/n", nndxin);    if (nndxin)    {        for (i = 0; i < nndxin; i++)        {            block2 = init_index(ndxinfiles[i], &gnames2);            srenew(gnames, block->nr+block2->nr);            for (j = 0; j < block2->nr; j++)            {                gnames[block->nr+j] = gnames2[j];            }            sfree(gnames2);            merge_blocks(block, block2);            sfree(block2->a);            sfree(block2->index);/*       done_block(block2); */            sfree(block2);        }    }    else    {        snew(gnames, 1);        analyse(atoms, block, &gnames, FALSE, TRUE);    }    if (!bNatoms)    {        natoms = block2natoms(block);        printf("Counted atom numbers up to %d in index file/n", natoms);    }    edit_index(natoms, atoms, x, block, &gnames, bVerbose);    write_index(ndxoutfile, block, gnames, bDuplicate, natoms);    return 0;}

static void dd_pmeredist_pos_coeffs(struct gmx_pme_t *pme,                                    int n, gmx_bool bX, rvec *x, real *data,                                    pme_atomcomm_t *atc){    int *commnode, *buf_index;    int  nnodes_comm, i, nsend, local_pos, buf_pos, node, scount, rcount;    commnode  = atc->node_dest;    buf_index = atc->buf_index;    nnodes_comm = min(2*atc->maxshift, atc->nslab-1);    nsend = 0;    for (i = 0; i < nnodes_comm; i++)    {        buf_index[commnode[i]] = nsend;        nsend                 += atc->count[commnode[i]];    }    if (bX)    {        if (atc->count[atc->nodeid] + nsend != n)        {            gmx_fatal(FARGS, "%d particles communicated to PME rank %d are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension %c./n"                      "This usually means that your system is not well equilibrated.",                      n - (atc->count[atc->nodeid] + nsend),                      pme->nodeid, 'x'+atc->dimind);        }        if (nsend > pme->buf_nalloc)        {            pme->buf_nalloc = over_alloc_dd(nsend);            srenew(pme->bufv, pme->buf_nalloc);            srenew(pme->bufr, pme->buf_nalloc);        }        atc->n = atc->count[atc->nodeid];        for (i = 0; i < nnodes_comm; i++)        {            scount = atc->count[commnode[i]];            /* Communicate the count */            if (debug)            {                fprintf(debug, "dimind %d PME rank %d send to rank %d: %d/n",                        atc->dimind, atc->nodeid, commnode[i], scount);            }            pme_dd_sendrecv(atc, FALSE, i,                            &scount, sizeof(int),                            &atc->rcount[i], sizeof(int));            atc->n += atc->rcount[i];        }        pme_realloc_atomcomm_things(atc);    }    local_pos = 0;    for (i = 0; i < n; i++)    {        node = atc->pd[i];        if (node == atc->nodeid)        {            /* Copy direct to the receive buffer */            if (bX)            {                copy_rvec(x[i], atc->x[local_pos]);            }            atc->coefficient[local_pos] = data[i];            local_pos++;        }        else        {            /* Copy to the send buffer */            if (bX)            {                copy_rvec(x[i], pme->bufv[buf_index[node]]);            }            pme->bufr[buf_index[node]] = data[i];            buf_index[node]++;        }    }    buf_pos = 0;    for (i = 0; i < nnodes_comm; i++)    {        scount = atc->count[commnode[i]];        rcount = atc->rcount[i];        if (scount > 0 || rcount > 0)        {            if (bX)            {                /* Communicate the coordinates */                pme_dd_sendrecv(atc, FALSE, i,                                pme->bufv[buf_pos], scount*sizeof(rvec),                                atc->x[local_pos], rcount*sizeof(rvec));            }            /* Communicate the coefficients */            pme_dd_sendrecv(atc, FALSE, i,                            pme->bufr+buf_pos, scount*sizeof(real),                            atc->coefficient+local_pos, rcount*sizeof(real));            buf_pos   += scount;            local_pos += atc->rcount[i];//.........这里部分代码省略.........

static int split_chain(t_atoms *atoms, rvec *x,                       int sel_nr, t_blocka *block, char ***gn){    char    buf[STRLEN];    int     j, nchain;    atom_id i, a, natoms, *start = NULL, *end = NULL, ca_start, ca_end;    rvec    vec;    natoms   = atoms->nr;    nchain   = 0;    ca_start = 0;    while (ca_start < natoms)    {        while ((ca_start < natoms) && strcmp(*atoms->atomname[ca_start], "CA"))        {            ca_start++;        }        if (ca_start < natoms)        {            srenew(start, nchain+1);            srenew(end, nchain+1);            start[nchain] = ca_start;            while ((start[nchain] > 0) &&                   (atoms->atom[start[nchain]-1].resind ==                    atoms->atom[ca_start].resind))            {                start[nchain]--;            }            i = ca_start;            do            {                ca_end = i;                do                {                    i++;                }                while ((i < natoms) && strcmp(*atoms->atomname[i], "CA"));                if (i < natoms)                {                    rvec_sub(x[ca_end], x[i], vec);                }            }            while ((i < natoms) && (norm(vec) < 0.45));            end[nchain] = ca_end;            while ((end[nchain]+1 < natoms) &&                   (atoms->atom[end[nchain]+1].resind == atoms->atom[ca_end].resind))            {                end[nchain]++;            }            ca_start = end[nchain]+1;            nchain++;        }    }    if (nchain == 1)    {        printf("Found 1 chain, will not split/n");    }    else    {        printf("Found %d chains/n", nchain);    }    for (j = 0; j < nchain; j++)    {        printf("%d:%6d atoms (%d to %d)/n",               j+1, end[j]-start[j]+1, start[j]+1, end[j]+1);    }    if (nchain > 1)    {        srenew(block->a, block->nra+block->index[sel_nr+1]-block->index[sel_nr]);        for (j = 0; j < nchain; j++)        {            block->nr++;            srenew(block->index, block->nr+1);            srenew(*gn, block->nr);            sprintf(buf, "%s_chain%d", (*gn)[sel_nr], j+1);            (*gn)[block->nr-1] = strdup(buf);            for (i = block->index[sel_nr]; i < block->index[sel_nr+1]; i++)            {                a = block->a[i];                if ((a >= start[j]) && (a <= end[j]))                {                    block->a[block->nra] = a;                    block->nra++;                }            }            block->index[block->nr] = block->nra;            if (block->index[block->nr-1] == block->index[block->nr])            {                remove_group(block->nr-1, NOTSET, block, gn);            }        }    }    sfree(start);    sfree(end);    return nchain;//.........这里部分代码省略.........

static void update_topol(const char *topinout, int p_num, int n_num,                         const char *p_name, const char *n_name, char *grpname){    FILE    *fpin, *fpout;    char     buf[STRLEN], buf2[STRLEN], *temp, **mol_line = NULL;    int      line, i, nmol_line, sol_line, nsol_last;    gmx_bool bMolecules;    char     temporary_filename[STRLEN];    printf("/nProcessing topology/n");    fpin  = gmx_ffopen(topinout, "r");    std::strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);    fpout = gmx_fopen_temporary(temporary_filename);    line       = 0;    bMolecules = FALSE;    nmol_line  = 0;    sol_line   = -1;    nsol_last  = -1;    while (fgets(buf, STRLEN, fpin))    {        line++;        std::strcpy(buf2, buf);        if ((temp = std::strchr(buf2, '/n')) != NULL)        {            temp[0] = '/0';        }        ltrim(buf2);        if (buf2[0] == '[')        {            buf2[0] = ' ';            if ((temp = std::strchr(buf2, '/n')) != NULL)            {                temp[0] = '/0';            }            rtrim(buf2);            if (buf2[std::strlen(buf2)-1] == ']')            {                buf2[std::strlen(buf2)-1] = '/0';                ltrim(buf2);                rtrim(buf2);                bMolecules = (gmx_strcasecmp(buf2, "molecules") == 0);            }            fprintf(fpout, "%s", buf);        }        else if (!bMolecules)        {            fprintf(fpout, "%s", buf);        }        else        {            /* Check if this is a line with solvent molecules */            sscanf(buf, "%s", buf2);            if (gmx_strcasecmp(buf2, grpname) == 0)            {                sol_line = nmol_line;                sscanf(buf, "%*s %d", &nsol_last);            }            /* Store this molecules section line */            srenew(mol_line, nmol_line+1);            mol_line[nmol_line] = gmx_strdup(buf);            nmol_line++;        }    }    gmx_ffclose(fpin);    if (sol_line == -1)    {        gmx_ffclose(fpout);        gmx_fatal(FARGS, "No line with moleculetype '%s' found the [ molecules ] section of file '%s'", grpname, topinout);    }    if (nsol_last < p_num+n_num)    {        gmx_ffclose(fpout);        gmx_fatal(FARGS, "The last entry for moleculetype '%s' in the [ molecules ] section of file '%s' has less solvent molecules (%d) than were replaced (%d)", grpname, topinout, nsol_last, p_num+n_num);    }    /* Print all the molecule entries */    for (i = 0; i < nmol_line; i++)    {        if (i != sol_line)        {            fprintf(fpout, "%s", mol_line[i]);        }        else        {            printf("Replacing %d solute molecules in topology file (%s) "                   " by %d %s and %d %s ions./n",                   p_num+n_num, topinout, p_num, p_name, n_num, n_name);            nsol_last -= p_num + n_num;            if (nsol_last > 0)            {                fprintf(fpout, "%-10s  %d/n", grpname, nsol_last);            }            if (p_num > 0)            {                fprintf(fpout, "%-15s  %d/n", p_name, p_num);            }            if (n_num > 0)            {//.........这里部分代码省略.........

/*! * /param[in,out] g      Index group structure. * * Resizes the memory allocated for holding the indices such that the * current contents fit. */voidgmx_ana_index_squeeze(gmx_ana_index_t *g){    srenew(g->index, g->isize);    g->nalloc_index = g->isize;}

//.........这里部分代码省略.........    {        rng = gmx_rng_init(gmx_rng_make_seed());    }    else    {        rng = gmx_rng_init(seed);    }    bTRX = ftp2bSet(efTRX, NFILE, fnm);    nx   = (int)(nrbox[XX]+0.5);    ny   = (int)(nrbox[YY]+0.5);    nz   = (int)(nrbox[ZZ]+0.5);    if ((nx <= 0) || (ny <= 0) || (nz <= 0))    {        gmx_fatal(FARGS, "Number of boxes (-nbox) should be larger than zero");    }    if ((nmolat <= 0) && bShuffle)    {        gmx_fatal(FARGS, "Can not shuffle if the molecules only have %d atoms",                  nmolat);    }    vol = nx*ny*nz; /* calculate volume in grid points (= nr. molecules) */    t_topology *top;    snew(top, 1);    atoms = &top->atoms;    read_tps_conf(opt2fn("-f", NFILE, fnm), top, &ePBC, &x, &v, box, FALSE);    natoms = atoms->nr;    nres   = atoms->nres;          /* nr of residues in one element? */    /* make space for all the atoms */    add_t_atoms(atoms, natoms*(vol-1), nres*(vol-1));    srenew(x, natoms*vol);         /* get space for coordinates of all atoms */    srenew(v, natoms*vol);         /* velocities. not really needed? */    snew(xrot, natoms);            /* get space for rotation matrix? */    snew(vrot, natoms);    if (bTRX)    {        if (!read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &xx, boxx))        {            gmx_fatal(FARGS, "No atoms in trajectory %s", ftp2fn(efTRX, NFILE, fnm));        }    }    else    {        snew(xx, natoms);        for (i = 0; i < natoms; i++)        {            copy_rvec(x[i], xx[i]);        }    }    for (k = 0; (k < nz); k++)     /* loop over all gridpositions    */    {        shift[ZZ] = k*(dist[ZZ]+box[ZZ][ZZ]);        for (j = 0; (j < ny); j++)        {            shift[YY] = j*(dist[YY]+box[YY][YY])+k*box[ZZ][YY];            for (i = 0; (i < nx); i++)            {                shift[XX] = i*(dist[XX]+box[XX][XX])+j*box[YY][XX]+k*box[ZZ][XX];

/*! * /param[in,out] t    Output block. * /param[in]  top  Topology structure *   (only used if /p type is /ref INDEX_RES or /ref INDEX_MOL, can be NULL *   otherwise). * /param[in]  g    Index group *   (can be NULL if /p type is /ref INDEX_UNKNOWN). * /param[in]  type Type of partitioning to make. * /param[in]  bComplete *   If true, the index group is expanded to include any residue/molecule *   (depending on /p type) that is partially contained in the group. *   If /p type is not INDEX_RES or INDEX_MOL, this has no effect. * * /p m should have been initialized somehow (calloc() is enough) unless * /p type is INDEX_UNKNOWN. * /p g should be sorted. */voidgmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,                         e_index_t type, bool bComplete){    int      i, j, ai;    int      id, cur;    if (type == INDEX_UNKNOWN)    {        t->nr           = 1;        snew(t->index, 2);        t->nalloc_index = 2;        t->index[0]     = 0;        t->index[1]     = 0;        t->nra          = 0;        t->a            = NULL;        t->nalloc_a     = 0;        return;    }    /* bComplete only does something for INDEX_RES or INDEX_MOL, so turn it     * off otherwise. */    if (type != INDEX_RES && type != INDEX_MOL)    {        bComplete = false;    }    /* Allocate memory for the atom array and fill it unless we are using     * completion. */    if (bComplete)    {        t->nra = 0;        /* We may allocate some extra memory here because we don't know in         * advance how much will be needed. */        if (t->nalloc_a < top->atoms.nr)        {            srenew(t->a, top->atoms.nr);            t->nalloc_a = top->atoms.nr;        }    }    else    {        t->nra      = g->isize;        if (t->nalloc_a < g->isize)        {            srenew(t->a, g->isize);            t->nalloc_a = g->isize;        }        std::memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));    }    /* Allocate memory for the block index. We don't know in advance     * how much will be needed, so we allocate some extra and free it in the     * end. */    if (t->nalloc_index < g->isize + 1)    {        srenew(t->index, g->isize + 1);        t->nalloc_index = g->isize + 1;    }    /* Clear counters */    t->nr = 0;    j     = 0; /* j is used by residue completion for the first atom not stored */    id    = cur = -1;    for (i = 0; i < g->isize; ++i)    {        ai = g->index[i];        /* Find the ID number of the atom/residue/molecule corresponding to         * atom ai. */        switch (type)        {            case INDEX_ATOM:                id = ai;                break;            case INDEX_RES:                id = top->atoms.atom[ai].resind;                break;            case INDEX_MOL:                while (ai >= top->mols.index[id+1])                {                    id++;                }                break;            case INDEX_UNKNOWN: /* Should not occur */            case INDEX_ALL://.........这里部分代码省略.........

static int low_fflib_search_file_end(const char *ffdir,                                     gmx_bool    bAddCWD,                                     const char *file_end,                                     gmx_bool    bFatalError,                                     char     ***filenames,                                     char     ***filenames_short){    char **fns = NULL, **fns_short = NULL;    int    n   = 0;    try    {        std::vector<std::string> libPaths;        bool                     bEnvIsSet = false;        if (ffdir != NULL)        {            /* Search ffdir in current dir and library dirs */            libPaths.push_back(gmxlibfn(ffdir));        }        else        {            /* GMXLIB can be a path now */            if (bAddCWD)            {                libPaths.push_back(".");            }            const char *lib = getenv("GMXLIB");            if (lib != NULL)            {                bEnvIsSet = true;                gmx::Path::splitPathEnvironment(lib, &libPaths);            }            else            {                libPaths.push_back(gmx::getProgramContext().defaultLibraryDataPath());            }        }        const int len_fe = strlen(file_end);        std::vector<std::string>::const_iterator i;        for (i = libPaths.begin(); i != libPaths.end(); ++i)        {            const char      *dir = i->c_str();            gmx_directory_t  dirhandle;            const int        rc  = gmx_directory_open(&dirhandle, dir);            if (rc == 0)            {                char nextname[STRLEN];                int  n_thisdir = 0;                while (gmx_directory_nextfile(dirhandle, nextname, STRLEN-1) == 0)                {                    nextname[STRLEN-1] = 0;                    if (debug)                    {                        fprintf(debug, "dir '%s' %d file '%s'/n",                                dir, n_thisdir, nextname);                    }                    const int len_name = strlen(nextname);                    /* What about case sensitivity? */                    if (len_name >= len_fe &&                        strcmp(nextname+len_name-len_fe, file_end) == 0)                    {                        char fn_dir[GMX_PATH_MAX];                        /* We have a match */                        srenew(fns, n+1);                        sprintf(fn_dir, "%s%c%s", dir, DIR_SEPARATOR, nextname);                        /* Copy the file name, possibly including the path. */                        fns[n] = gmx_strdup(fn_dir);                        if (ffdir == NULL)                        {                            /* We are searching in a path.                             * Use the relative path when we use share/top                             * from the installation.                             * Add the full path when we use the current                             * working directory of GMXLIB.                             */                            srenew(fns_short, n+1);                            if (strcmp(dir, ".") == 0 || bEnvIsSet)                            {                                fns_short[n] = gmx_strdup(fn_dir);                            }                            else                            {                                fns_short[n] = gmx_strdup(nextname);                            }                        }                        n++;                        n_thisdir++;                    }                }                gmx_directory_close(dirhandle);                sort_filenames(n_thisdir,                               fns+n-n_thisdir,                               fns_short == NULL ? NULL : fns_short+n-n_thisdir);            }        }//.........这里部分代码省略.........

int sscan_list(int *list[], const char *str, const char *listname){    /*this routine scans a string of the form 1,3-6,9 and returns the       selected numbers (in this case 1 3 4 5 6 9) in NULL-terminated array of integers.       memory for this list will be allocated  in this routine -- sscan_list expects *list to       be a NULL-Pointer       listname is a string used in the errormessage*/    int   i, istep;    char  c;    char *pos, *startpos, *step;    int   n = strlen(str);    /*enums to define the different lexical stati */    enum {        sBefore, sNumber, sMinus, sRange, sZero, sSmaller, sError, sSteppedRange    };    int   status     = sBefore; /*status of the deterministic automat to scan str   */    int   number     = 0;    int   end_number = 0;    char *start = NULL; /*holds the string of the number behind a ','*/    char *end   = NULL; /*holds the string of the number behind a '-' */    int   nvecs = 0;    /* counts the number of vectors in the list*/    step = NULL;    snew(pos, n+4);    startpos = pos;    strcpy(pos, str);    pos[n]   = ',';    pos[n+1] = '1';    pos[n+2] = '/0';    *list = NULL;    while ((c = *pos) != 0)    {        switch (status)        {            /* expect a number */            case sBefore: if (isdigit(c))                {                    start  = pos;                    status = sNumber;                    break;                }                else                {                    status = sError;                } break;            /* have read a number, expect ',' or '-' */            case sNumber: if (c == ',')                {                    /*store number*/                    srenew(*list, nvecs+1);                    (*list)[nvecs++] = number = strtol(start, NULL, 10);                    status           = sBefore;                    if (number == 0)                    {                        status = sZero;                    }                    break;                }                else if (c == '-')                {                    status = sMinus; break;                }                else if (isdigit(c))                {                    break;                }                else                {                    status = sError;                } break;            /* have read a '-' -> expect a number */            case sMinus:                if (isdigit(c))                {                    end    = pos;                    status = sRange; break;                }                else                {                    status = sError;                } break;            case sSteppedRange:                if (isdigit(c))                {                    if (step)                    {                        status = sError; break;                    }//.........这里部分代码省略.........