static t_bbb *mk_bonds(int natoms,rvec x[],real odist,
		       gmx_bool bPBC,matrix box)
{
  real  od2 = odist*odist+1e-5;
  t_pbc pbc;
  t_bbb *bbb;
  int   i,j;
  rvec  dx;
  
  if (bPBC)
    set_pbc(&pbc,box);
  snew(bbb,natoms);
  for(i=0; (i<natoms); i++) {
    for(j=i+1; (j<natoms); j++) {
      if (bPBC)
	pbc_dx(&pbc,x[i],x[j],dx);
      else
	rvec_sub(x[i],x[j],dx);
      if (iprod(dx,dx) <= od2) {
	bbb[i].aa[bbb[i].n++] = j;
	bbb[j].aa[bbb[j].n++] = i;
      }
    }
  }
  if (debug) 
#define PRB(nn) (bbb[(i)].n >= (nn)) ? bbb[i].aa[nn-1] : -1
    for(i=0; (i<natoms); i++)
      fprintf(debug,"bonds from %3d:  %d %d %d %d\n",
	      i,PRB(1),PRB(2),PRB(3),PRB(4));
#undef PRB
  return bbb;
}

/* This is a (maybe) slow workaround to avoid the neighbor searching in addconf.c, which
 * leaks memory (May 2012). The function could be deleted as soon as the momory leaks
 * in addconf.c are fixed.
 * However, when inserting a small molecule in a system containing not too many atoms,
 * allPairsDistOk is probably even faster than addconf.c
 */
static gmx_bool
allPairsDistOk(t_atoms *atoms, rvec *x, real *r,
               int ePBC, matrix box,
               t_atoms *atoms_insrt, rvec *x_n, real *r_insrt)
{
    int   i, j;
    rvec  dx;
    real  n2, r2;
    t_pbc pbc;

    set_pbc(&pbc, ePBC, box);
    for (i = 0; i < atoms->nr; i++)
    {
        for (j = 0; j < atoms_insrt->nr; j++)
        {
            pbc_dx(&pbc, x[i], x_n[j], dx);
            n2 = norm2(dx);
            r2 = sqr(r[i]+r_insrt[j]);
            if (n2 < r2)
            {
                return FALSE;
            }
        }
    }
    return TRUE;
}

void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,
                   const output_env_t oenv)
{
    pull_params_t *pull;
    struct pull_t *pull_work;
    t_mdatoms     *md;
    t_pbc          pbc;
    int            c;
    double         t_start;

    pull      = ir->pull;
    pull_work = init_pull(NULL, pull, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
    md        = init_mdatoms(NULL, mtop, ir->efep);
    atoms2md(mtop, ir, 0, NULL, mtop->natoms, md);
    if (ir->efep)
    {
        update_mdatoms(md, lambda);
    }

    set_pbc(&pbc, ir->ePBC, box);

    t_start = ir->init_t + ir->init_step*ir->delta_t;

    pull_calc_coms(NULL, pull_work, md, &pbc, t_start, x, NULL);

    fprintf(stderr, "Pull group  natoms  pbc atom  distance at start  reference at t=0\n");
    for (c = 0; c < pull->ncoord; c++)
    {
        t_pull_coord *pcrd;
        t_pull_group *pgrp0, *pgrp1;
        double        value;
        real          init = 0;

        pcrd  = &pull->coord[c];

        pgrp0 = &pull->group[pcrd->group[0]];
        pgrp1 = &pull->group[pcrd->group[1]];
        fprintf(stderr, "%8d  %8d  %8d\n",
                pcrd->group[0], pgrp0->nat, pgrp0->pbcatom+1);
        fprintf(stderr, "%8d  %8d  %8d ",
                pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1);

        if (pcrd->bStart)
        {
            init       = pcrd->init;
            pcrd->init = 0;
        }

        get_pull_coord_value(pull_work, c, &pbc, &value);
        fprintf(stderr, " %10.3f nm", value);

        if (pcrd->bStart)
        {
            pcrd->init = value + init;
        }
        fprintf(stderr, "     %10.3f nm\n", pcrd->init);
    }

    finish_pull(pull_work);
}

static void calc_dist_tot(int nind, atom_id index[], rvec x[],
                          int ePBC, matrix box,
                          real **d, real **dtot, real **dtot2,
                          gmx_bool bNMR, real **dtot1_3, real **dtot1_6)
{
    int      i, j;
    real    *xi;
    real     temp, temp2, temp1_3;
    rvec     dx;
    t_pbc    pbc;

    set_pbc(&pbc, ePBC, box);
    for (i = 0; (i < nind-1); i++)
    {
        xi = x[index[i]];
        for (j = i+1; (j < nind); j++)
        {
            pbc_dx(&pbc, xi, x[index[j]], dx);
            temp2        = norm2(dx);
            temp         = std::sqrt(temp2);
            d[i][j]      = temp;
            dtot[i][j]  += temp;
            dtot2[i][j] += temp2;
            if (bNMR)
            {
                temp1_3        = 1.0/(temp*temp2);
                dtot1_3[i][j] += temp1_3;
                dtot1_6[i][j] += temp1_3*temp1_3;
            }
        }
    }
}

void mk_bonds(int nnm, t_nm2type nmt[],
              t_atoms *atoms, rvec x[], t_params *bond, int nbond[],
              gmx_bool bPBC, matrix box)
{
    t_param b;
    int     i, j;
    t_pbc   pbc;
    rvec    dx;
    real    dx2;

    for (i = 0; (i < MAXATOMLIST); i++)
    {
        b.a[i] = -1;
    }
    for (i = 0; (i < MAXFORCEPARAM); i++)
    {
        b.c[i] = 0.0;
    }

    if (bPBC)
    {
        set_pbc(&pbc, -1, box);
    }
    for (i = 0; (i < atoms->nr); i++)
    {
        if ((i % 10) == 0)
        {
            fprintf(stderr, "\ratom %d", i);
        }
        for (j = i+1; (j < atoms->nr); j++)
        {
            if (bPBC)
            {
                pbc_dx(&pbc, x[i], x[j], dx);
            }
            else
            {
                rvec_sub(x[i], x[j], dx);
            }

            dx2 = iprod(dx, dx);
            if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
                        sqrt(dx2)))
            {
                b.AI = i;
                b.AJ = j;
                b.C0 = sqrt(dx2);
                add_param_to_list (bond, &b);
                nbond[i]++;
                nbond[j]++;
                if (debug)
                {
                    fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
                            *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
                }
            }
        }
    }
    fprintf(stderr, "\ratom %d\n", i);
}

static void chk_bonds(t_idef *idef, int ePBC, rvec *x, matrix box, real tol)
{
    int   ftype, i, k, ai, aj, type;
    real  b0, blen, deviation, devtot;
    t_pbc pbc;
    rvec  dx;

    devtot = 0;
    set_pbc(&pbc, ePBC, box);
    for (ftype = 0; (ftype < F_NRE); ftype++)
    {
        if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
        {
            for (k = 0; (k < idef->il[ftype].nr); )
            {
                type = idef->il[ftype].iatoms[k++];
                ai   = idef->il[ftype].iatoms[k++];
                aj   = idef->il[ftype].iatoms[k++];
                b0   = 0;
                switch (ftype)
                {
                    case F_BONDS:
                        b0 = idef->iparams[type].harmonic.rA;
                        break;
                    case F_G96BONDS:
                        b0 = sqrt(idef->iparams[type].harmonic.rA);
                        break;
                    case F_MORSE:
                        b0 = idef->iparams[type].morse.b0A;
                        break;
                    case F_CUBICBONDS:
                        b0 = idef->iparams[type].cubic.b0;
                        break;
                    case F_CONSTR:
                        b0 = idef->iparams[type].constr.dA;
                        break;
                    default:
                        break;
                }
                if (b0 != 0)
                {
                    pbc_dx(&pbc, x[ai], x[aj], dx);
                    blen      = norm(dx);
                    deviation = sqr(blen-b0);
                    if (sqrt(deviation/sqr(b0) > tol))
                    {
                        fprintf(stderr, "Distance between atoms %d and %d is %.3f, should be %.3f\n", ai+1, aj+1, blen, b0);
                    }
                }
            }
        }
    }
}

static void calc_mj(t_topology top, int ePBC, matrix box, gmx_bool bNoJump, int isize, int index0[], \
                    rvec fr[], rvec mj, real mass2[], real qmol[])
{

    int   i, j, k, l;
    rvec  tmp;
    rvec  center;
    rvec  mt1, mt2;
    t_pbc pbc;


    calc_box_center(ecenterRECT, box, center);

    if (!bNoJump)
    {
        set_pbc(&pbc, ePBC, box);
    }

    clear_rvec(tmp);


    for (i = 0; i < isize; i++)
    {
        clear_rvec(mt1);
        clear_rvec(mt2);
        k = top.mols.index[index0[i]];
        l = top.mols.index[index0[i+1]];
        for (j = k; j < l; j++)
        {
            svmul(mass2[j], fr[j], tmp);
            rvec_inc(mt1, tmp);
        }

        if (bNoJump)
        {
            svmul(qmol[k], mt1, mt1);
        }
        else
        {
            pbc_dx(&pbc, mt1, center, mt2);
            svmul(qmol[k], mt2, mt1);
        }

        rvec_inc(mj, mt1);

    }

}

void setStateDependentAwhParams(AwhParams *awhParams,
                                const pull_params_t *pull_params, pull_t *pull_work,
                                const matrix box,  int ePBC,
                                const t_grpopts *inputrecGroupOptions, warninp_t wi)
{
    /* The temperature is not really state depenendent but is not known
     * when read_awhParams is called (in get ir).
     * It is known first after do_index has been called in grompp.cpp.
     */
    if (inputrecGroupOptions->ref_t == NULL ||
        inputrecGroupOptions->ref_t[0] <= 0)
    {
        gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
    }
    for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
    {
        if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
        {
            gmx_fatal(FARGS, "AWH biasing is currently only supported for identical temperatures for all temperature coupling groups");
        }
    }

    t_pbc          pbc;
    set_pbc(&pbc, ePBC, box);

    for (int k = 0; k < awhParams->numBias; k++)
    {
        AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
        for (int d = 0; d < awhBiasParams->ndim; d++)
        {
            AwhDimParams *dimParams = &awhBiasParams->dimParams[d];

            /* The periodiciy of the AWH grid in certain cases depends on the simulation box */
            dimParams->period = get_pull_coord_period(pull_params, dimParams->coordIndex, box);

            /* The initial coordinate value, converted to external user units. */
            dimParams->coordValueInit =
                get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);

            t_pull_coord *pullCoord = &pull_params->coord[dimParams->coordIndex];
            dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(pullCoord);
        }
    }
    checkInputConsistencyInterval(awhParams, wi);

    /* Register AWH as external potential with pull to check consistency. */
    Awh::registerAwhWithPull(*awhParams, pull_work);
}

void calc_angles_dihs(t_params *ang, t_params *dih, rvec x[], gmx_bool bPBC,
                      matrix box)
{
    int    i, ai, aj, ak, al, t1, t2, t3;
    rvec   r_ij, r_kj, r_kl, m, n;
    real   sign, th, costh, ph;
    t_pbc  pbc;

    if (bPBC)
    {
        set_pbc(&pbc, epbcXYZ, box);
    }
    if (debug)
    {
        pr_rvecs(debug, 0, "X2TOP", box, DIM);
    }
    for (i = 0; (i < ang->nr); i++)
    {
        ai = ang->param[i].AI;
        aj = ang->param[i].AJ;
        ak = ang->param[i].AK;
        th = RAD2DEG*bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : NULL,
                                r_ij, r_kj, &costh, &t1, &t2);
        if (debug)
        {
            fprintf(debug, "X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
                    ai, aj, ak, norm(r_ij), norm(r_kj), th);
        }
        ang->param[i].C0 = th;
    }
    for (i = 0; (i < dih->nr); i++)
    {
        ai = dih->param[i].AI;
        aj = dih->param[i].AJ;
        ak = dih->param[i].AK;
        al = dih->param[i].AL;
        ph = RAD2DEG*dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : NULL,
                               r_ij, r_kj, r_kl, m, n, &sign, &t1, &t2, &t3);
        if (debug)
        {
            fprintf(debug, "X2TOP: ai=%3d aj=%3d ak=%3d al=%3d r_ij=%8.3f r_kj=%8.3f r_kl=%8.3f ph=%8.3f\n",
                    ai, aj, ak, al, norm(r_ij), norm(r_kj), norm(r_kl), ph);
        }
        dih->param[i].C0 = ph;
    }
}

static void calc_mat(int nres, int natoms, int rndx[],
                     rvec x[], atom_id *index,
                     real trunc, real **mdmat, int **nmat, int ePBC, matrix box)
{
    int   i, j, resi, resj;
    real  trunc2, r, r2;
    t_pbc pbc;
    rvec  ddx;

    set_pbc(&pbc, ePBC, box);
    trunc2 = sqr(trunc);
    for (resi = 0; (resi < nres); resi++)
    {
        for (resj = 0; (resj < nres); resj++)
        {
            mdmat[resi][resj] = FARAWAY;
        }
    }
    for (i = 0; (i < natoms); i++)
    {
        resi = rndx[i];
        for (j = i+1; (j < natoms); j++)
        {
            resj = rndx[j];
            pbc_dx(&pbc, x[index[i]], x[index[j]], ddx);
            r2 = norm2(ddx);
            if (r2 < trunc2)
            {
                nmat[resi][j]++;
                nmat[resj][i]++;
            }
            mdmat[resi][resj] = min(r2, mdmat[resi][resj]);
        }
    }

    for (resi = 0; (resi < nres); resi++)
    {
        mdmat[resi][resi] = 0;
        for (resj = resi+1; (resj < nres); resj++)
        {
            r                 = sqrt(mdmat[resi][resj]);
            mdmat[resi][resj] = r;
            mdmat[resj][resi] = r;
        }
    }
}

const char *
put_atoms_in_compact_unitcell(int ePBC,int ecenter,matrix box,
                              int natoms,rvec x[])
                              {
    t_pbc pbc;
    rvec box_center,dx;
    int  i;

    set_pbc(&pbc,ePBC,box);
    calc_box_center(ecenter,box,box_center);
    for(i=0; i<natoms; i++) {
        pbc_dx(&pbc,x[i],box_center,dx);
        rvec_add(box_center,dx,x[i]);
    }

    return pbc.bLimitDistance ?
        "WARNING: Could not put all atoms in the compact unitcell\n"
        : NULL;
                              }

static void calc_dist(int nind, int index[], const rvec x[], int ePBC, matrix box,
                      real **d)
{
    int      i, j;
    rvec     dx;
    real     temp2;
    t_pbc    pbc;

    set_pbc(&pbc, ePBC, box);
    for (i = 0; (i < nind-1); i++)
    {
        const real *xi = x[index[i]];
        for (j = i+1; (j < nind); j++)
        {
            pbc_dx(&pbc, xi, x[index[j]], dx);
            temp2   = norm2(dx);
            d[i][j] = std::sqrt(temp2);
        }
    }
}

void
partition_absorbing_box( grid_t * g,
                         double gx0, double gy0, double gz0,
                         double gx1, double gy1, double gz1,
                         int gnx, int gny, int gnz,
                         int gpx, int gpy, int gpz,
                         int pbc ) {
  int px, py, pz; 

  partition_periodic_box( g,
                          gx0, gy0, gz0,
                          gx1, gy1, gz1,
                          gnx, gny, gnz,
                          gpx, gpy, gpz );

  // Override periodic boundary conditions

  RANK_TO_INDEX( world_rank, px,py,pz );

  if( px==0 && gnx>1 ) { 
    set_fbc(g,BOUNDARY(-1,0,0),absorb_fields);
    set_pbc(g,BOUNDARY(-1,0,0),pbc);
  } 

  if( px==gpx-1 && gnx>1 ) {
    set_fbc(g,BOUNDARY( 1,0,0),absorb_fields);
    set_pbc(g,BOUNDARY( 1,0,0),pbc);
  }

  if( py==0 && gny>1 ) { 
    set_fbc(g,BOUNDARY(0,-1,0),absorb_fields);
    set_pbc(g,BOUNDARY(0,-1,0),pbc);
  } 

  if( py==gpy-1 && gny>1 ) {
    set_fbc(g,BOUNDARY(0, 1,0),absorb_fields);
    set_pbc(g,BOUNDARY(0, 1,0),pbc);
  }

  if( pz==0 && gnz>1 ) { 
    set_fbc(g,BOUNDARY(0,0,-1),absorb_fields);
    set_pbc(g,BOUNDARY(0,0,-1),pbc);
  } 

  if( pz==gpz-1 && gnz>1 ) {
    set_fbc(g,BOUNDARY(0,0, 1),absorb_fields);
    set_pbc(g,BOUNDARY(0,0, 1),pbc);
  }
}

void put_atoms_in_compact_unitcell(int ePBC, int ecenter, const matrix box,
                                   int natoms, rvec x[])
{
    t_pbc pbc;
    rvec  box_center, dx;
    int   i;

    set_pbc(&pbc, ePBC, box);

    if (pbc.ePBCDX == epbcdxUNSUPPORTED)
    {
        gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
    }

    calc_box_center(ecenter, box, box_center);
    for (i = 0; i < natoms; i++)
    {
        pbc_dx(&pbc, x[i], box_center, dx);
        rvec_add(box_center, dx, x[i]);
    }
}

//.........这里部分代码省略.........
      if (index[g][i]>max)
	max=index[g][i];
      if (index[g][i] >= top->atoms.nr)
	gmx_fatal(FARGS,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);
      mass[g]+=top->atoms.atom[index[g][i]].m;
    }
  }

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

  if (max>=natoms)
    gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);

  if (!bCutoff) {
    /* open output file */
    fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),
		  "Distance","Time (ps)","Distance (nm)",oenv);
    xvgr_legend(fp,4,leg,oenv);
  } else {
    ngrps = 1;
    if (bLifeTime)
      snew(contact_time,isize[1]);
  }
  if (ePBC != epbcNONE)
    snew(pbc,1);
  else
    pbc = NULL;
    
  gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
  do {
    /* initialisation for correct distance calculations */
    if (pbc) {
      set_pbc(pbc,ePBC,box);
      /* make molecules whole again */
      gmx_rmpbc(gpbc,natoms,box,x);
    }
    /* calculate center of masses */
    for(g=0;(g<ngrps);g++) {
      if (isize[g] == 1) {
	copy_rvec(x[index[g][0]],com[g]);
      } else {
	for(d=0;(d<DIM);d++) {
	  com[g][d]=0;
	  for(i=0;(i<isize[g]);i++) {
	    com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
	  }
	  com[g][d] /= mass[g];
	}
      }
    }
    
    if (!bCutoff) {
      /* write to output */
      fprintf(fp,"%12.7f ",t);
      for(g=0;(g<ngrps/2);g++) {
	if (pbc)
	  pbc_dx(pbc,com[2*g],com[2*g+1],dx);
	else
	  rvec_sub(com[2*g],com[2*g+1],dx);
	
	fprintf(fp,"%12.7f %12.7f %12.7f %12.7f",
		norm(dx),dx[XX],dx[YY],dx[ZZ]);
      }
      fprintf(fp,"\n");
    } else {

//.........这里部分代码省略.........
    for(i=0; i<n1; i++)
        snew(grid[i],n2);

    box1 = 0;
    box2 = 0;
    nfr = 0;
    do {
        if (!bRadial) {
            box1 += box[c1][c1];
            box2 += box[c2][c2];
            invcellvol = n1*n2;
            if (nmpower == -3)
                invcellvol /= det(box);
            else if (nmpower == -2)
                invcellvol /= box[c1][c1]*box[c2][c2];
            for(i=0; i<nindex; i++) {
                j = index[i];
                if ((!bXmin || x[j][cav] >= xmin) &&
                        (!bXmax || x[j][cav] <= xmax)) {
                    m1 = x[j][c1]/box[c1][c1];
                    if (m1 >= 1)
                        m1 -= 1;
                    if (m1 < 0)
                        m1 += 1;
                    m2 = x[j][c2]/box[c2][c2];
                    if (m2 >= 1)
                        m2 -= 1;
                    if (m2 < 0)
                        m2 += 1;
                    grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
                }
            }
        } else {
            set_pbc(&pbc,ePBC,box);
            for(i=0; i<2; i++) {
                if (gnx[i] == 1) {
                    /* One atom, just copy the coordinates */
                    copy_rvec(x[ind[i][0]],xcom[i]);
                } else {
                    /* Calculate the center of mass */
                    clear_rvec(xcom[i]);
                    mtot = 0;
                    for(j=0; j<gnx[i]; j++) {
                        k = ind[i][j];
                        m = top.atoms.atom[k].m;
                        for(l=0; l<DIM; l++)
                            xcom[i][l] += m*x[k][l];
                        mtot += m;
                    }
                    svmul(1/mtot,xcom[i],xcom[i]);
                }
            }
            pbc_dx(&pbc,xcom[1],xcom[0],direction);
            for(i=0; i<DIM; i++)
                center[i] = xcom[0][i] + 0.5*direction[i];
            unitv(direction,direction);
            for(i=0; i<nindex; i++) {
                j = index[i];
                pbc_dx(&pbc,x[j],center,dx);
                axial = iprod(dx,direction);
                r = sqrt(norm2(dx) - axial*axial);
                if (axial>=-amax && axial<amax && r<rmax) {
                    if (bMirror)
                        r += rmax;
                    grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
                }

void read_ang_dih(const char *trj_fn,
                  gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
                  int maxangstat, int angstat[],
                  int *nframes, real **time,
                  int isize, atom_id index[],
                  real **trans_frac,
                  real **aver_angle,
                  real *dih[],
                  const output_env_t oenv)
{
    t_pbc       *pbc;
    t_trxstatus *status;
    int          i, angind, natoms, total, teller;
    int          nangles, n_alloc;
    real         t, fraction, pifac, aa, angle;
    real        *angles[2];
    matrix       box;
    rvec        *x;
    int          cur = 0;
#define prev (1-cur)

    snew(pbc, 1);
    natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);

    if (bAngles)
    {
        nangles = isize/3;
        pifac   = M_PI;
    }
    else
    {
        nangles = isize/4;
        pifac   = 2.0*M_PI;
    }
    snew(angles[cur], nangles);
    snew(angles[prev], nangles);

    /* Start the loop over frames */
    total       = 0;
    teller      = 0;
    n_alloc     = 0;
    *time       = NULL;
    *trans_frac = NULL;
    *aver_angle = NULL;

    do
    {
        if (teller >= n_alloc)
        {
            n_alloc += 100;
            if (bSaveAll)
            {
                for (i = 0; (i < nangles); i++)
                {
                    srenew(dih[i], n_alloc);
                }
            }
            srenew(*time, n_alloc);
            srenew(*trans_frac, n_alloc);
            srenew(*aver_angle, n_alloc);
        }

        (*time)[teller] = t;

        if (pbc)
        {
            set_pbc(pbc, -1, box);
        }

        if (bAngles)
        {
            calc_angles(pbc, isize, index, angles[cur], x);
        }
        else
        {
            calc_dihs(pbc, isize, index, angles[cur], x);

            /* Trans fraction */
            fraction              = calc_fraction(angles[cur], nangles);
            (*trans_frac)[teller] = fraction;

            /* Change Ryckaert-Bellemans dihedrals to polymer convention
             * Modified 990913 by Erik:
             * We actually shouldn't change the convention, since it's
             * calculated from polymer above, but we change the intervall
             * from [-180,180] to [0,360].
             */
            if (bRb)
            {
                for (i = 0; (i < nangles); i++)
                {
                    if (angles[cur][i] <= 0.0)
                    {
                        angles[cur][i] += 2*M_PI;
                    }
                }
            }

            /* Periodicity in dihedral space... */
            if (bPBC)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                snew(khist, histbins);
            }

            do
            {
                clear_rvec(donvec);
                clear_rvec(accvec);
                clear_rvec(donpos);
                clear_rvec(accpos);
                for (i = 0; i < ndon / 2; i++)
                {
                    rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
                    rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
                    rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
                    rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
                }

                for (i = 0; i < nacc / 2; i++)
                {
                    rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
                    rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
                    rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
                    rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
                }

                unitv(donvec, donvec);
                unitv(accvec, accvec);

                svmul(1.0 / ndon, donpos, donpos);
                svmul(1.0 / nacc, accpos, accpos);

                if (bPBCdist)
                {
                    set_pbc(pbc, ePBC, fr.box);
                    pbc_dx(pbc, donpos, accpos, dist);
                }
                else
                {
                    rvec_sub(donpos, accpos, dist);
                }

                unitv(dist, distnorm);
                R       = norm(dist);
                kappa2  = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
                kappa2 *= kappa2;
                if (R0 > 0)
                {
                    Rfrac     = R/R0;
                    insteff   = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
                    insteffs += insteff;

                    if (bInstEffout)
                    {
                        fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
                    }
                }


                Rs      += R;
                kappa2s += kappa2;
                rkcount++;

                if (bRKout)
                {
                    fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
                }

void orient(int natom,rvec *x,rvec *v, rvec angle,matrix box)
{
  real longest,rij,rzi;
  int  i,j,m,max_i=0,max_j=0;
  rvec origin;
  int  temp;
  real alfa=0,beta=0,gamma=0;
  t_pbc pbc;

  set_pbc(&pbc,-1,box);

  /*first i am going to look for the longest atom-atom distance*/
  longest=dist2(&pbc,x[0],x[1]);
  i=0;
  j=1;
  for (i=0;(i<natom);i++) {
    for (j=0;(j<natom);j++) {
      rij=dist2(&pbc,x[i],x[j]);
      if (rij>longest) {
	max_i=i;
	max_j=j;
	longest=rij;
      }
    }
  }
  /* first check if x[max_i]<x[max_j] else swap*/
  if (x[max_i][2]>x[max_j][2]) {
    temp=max_i;
    max_i=max_j;
    max_j=temp;
  }
  
  /*set the origin to x[i]*/
  for(m=0;(m<DIM);m++) 
    origin[m]=x[max_i][m];
  for(i=0;(i<natom);i++)
    for(m=0;(m<DIM);m++)
      x[i][m]-=origin[m];
      
  /* calculate the rotation angles alfa(x_axis) and beta(y_axis)
   * the rotation angles must be calculated clockwise looking 
   * along the rotation axis to the origin*
   * alfa (x-axis)
   */
  alfa=atan(x[max_j][ZZ]/x[max_j][YY])-M_PI_2;
  beta=M_PI_2-atan(x[max_j][ZZ]/x[max_j][XX]);
  rotate_conf(natom,x,v,alfa,beta,gamma);
  
  /* now search the longest distance for rotation along the z_axis */
  longest=distance_to_z(x[0]);
  max_i=0;
  for (i=1;(i<natom);i++) {
    rzi=distance_to_z(x[i]);
    if (rzi>longest) {
      longest = rzi;
      max_i=i;
    }
  }
  gamma=atan(x[max_i][YY]/x[max_i][XX])-M_PI_2;
  rotate_conf(natom,x,v,0,0,gamma);
  angle[0]=alfa;
  angle[1]=beta;
  angle[2]=gamma;
} /*orient()*/

void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                       t_idef     *idef,    t_commrec  *cr,
                       t_nrnb     *nrnb,    gmx_wallcycle_t wcycle,
                       t_mdatoms  *md,
                       rvec       x[],      history_t  *hist,
                       rvec       f[],
                       rvec       f_longrange[],
                       gmx_enerdata_t *enerd,
                       t_fcdata   *fcd,
                       gmx_localtop_t *top,
                       gmx_genborn_t *born,
                       gmx_bool       bBornRadii,
                       matrix     box,
                       t_lambda   *fepvals,
                       real       *lambda,
                       t_graph    *graph,
                       t_blocka   *excl,
                       rvec       mu_tot[],
                       int        flags,
                       float      *cycles_pme)
{
    int         i, j;
    int         donb_flags;
    gmx_bool    bSB;
    int         pme_flags;
    matrix      boxs;
    rvec        box_size;
    t_pbc       pbc;
    real        dvdl_dum[efptNR], dvdl_nb[efptNR];

#ifdef GMX_MPI
    double  t0 = 0.0, t1, t2, t3; /* time measurement for coarse load balancing */
#endif

    set_pbc(&pbc, fr->ePBC, box);

    /* reset free energy components */
    for (i = 0; i < efptNR; i++)
    {
        dvdl_nb[i]  = 0;
        dvdl_dum[i] = 0;
    }

    /* Reset box */
    for (i = 0; (i < DIM); i++)
    {
        box_size[i] = box[i][i];
    }

    debug_gmx();

    /* do QMMM first if requested */
    if (fr->bQMMM)
    {
        enerd->term[F_EQM] = calculate_QMMM(cr, x, f, fr);
    }

    /* Call the short range functions all in one go. */

#ifdef GMX_MPI
    /*#define TAKETIME ((cr->npmenodes) && (fr->timesteps < 12))*/
#define TAKETIME FALSE
    if (TAKETIME)
    {
        MPI_Barrier(cr->mpi_comm_mygroup);
        t0 = MPI_Wtime();
    }
#endif

    if (ir->nwall)
    {
        /* foreign lambda component for walls */
        real dvdl_walls = do_walls(ir, fr, box, md, x, f, lambda[efptVDW],
                                   enerd->grpp.ener[egLJSR], nrnb);
        enerd->dvdl_lin[efptVDW] += dvdl_walls;
    }

    /* If doing GB, reset dvda and calculate the Born radii */
    if (ir->implicit_solvent)
    {
        wallcycle_sub_start(wcycle, ewcsNONBONDED);

        for (i = 0; i < born->nr; i++)
        {
            fr->dvda[i] = 0;
        }

        if (bBornRadii)
        {
            calc_gb_rad(cr, fr, ir, top, x, &(fr->gblist), born, md, nrnb);
        }

        wallcycle_sub_stop(wcycle, ewcsNONBONDED);
    }

    where();
    /* We only do non-bonded calculation with group scheme here, the verlet
     * calls are done from do_force_cutsVERLET(). */
    if (fr->cutoff_scheme == ecutsGROUP && (flags & GMX_FORCE_NONBONDED))
    {
//.........这里部分代码省略.........