C++ copy_rvec函数代码示例

static void make_cyl_refgrps(t_commrec *cr, struct pull_t *pull, t_mdatoms *md,
                             t_pbc *pbc, double t, rvec *x)
{
    /* The size and stride per coord for the reduction buffer */
    const int       stride = 9;
    int             c, i, ii, m, start, end;
    rvec            g_x, dx, dir;
    double          inv_cyl_r2;
    pull_comm_t    *comm;
    gmx_ga2la_t    *ga2la = NULL;

    comm = &pull->comm;

    if (comm->dbuf_cyl == NULL)
    {
        snew(comm->dbuf_cyl, pull->ncoord*stride);
    }

    if (cr && DOMAINDECOMP(cr))
    {
        ga2la = cr->dd->ga2la;
    }

    start = 0;
    end   = md->homenr;

    inv_cyl_r2 = 1.0/gmx::square(pull->params.cylinder_r);

    /* loop over all groups to make a reference group for each*/
    for (c = 0; c < pull->ncoord; c++)
    {
        pull_coord_work_t *pcrd;
        double             sum_a, wmass, wwmass;
        dvec               radf_fac0, radf_fac1;

        pcrd   = &pull->coord[c];

        sum_a  = 0;
        wmass  = 0;
        wwmass = 0;
        clear_dvec(radf_fac0);
        clear_dvec(radf_fac1);

        if (pcrd->params.eGeom == epullgCYL)
        {
            pull_group_work_t *pref, *pgrp, *pdyna;

            /* pref will be the same group for all pull coordinates */
            pref  = &pull->group[pcrd->params.group[0]];
            pgrp  = &pull->group[pcrd->params.group[1]];
            pdyna = &pull->dyna[c];
            copy_rvec(pcrd->vec, dir);
            pdyna->nat_loc = 0;

            /* We calculate distances with respect to the reference location
             * of this cylinder group (g_x), which we already have now since
             * we reduced the other group COM over the ranks. This resolves
             * any PBC issues and we don't need to use a PBC-atom here.
             */
            if (pcrd->params.rate != 0)
            {
                /* With rate=0, value_ref is set initially */
                pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
            }
            for (m = 0; m < DIM; m++)
            {
                g_x[m] = pgrp->x[m] - pcrd->vec[m]*pcrd->value_ref;
            }

            /* loop over all atoms in the main ref group */
            for (i = 0; i < pref->params.nat; i++)
            {
                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;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    atoms2md(top_global, inputrec, 0, NULL, top_global->natoms, mdatoms);
    update_mdatoms(mdatoms, inputrec->fepvals->init_lambda);

    snew(enerd, 1);
    init_enerdata(groups->grps[egcENER].nr, inputrec->fepvals->n_lambda, enerd);
    snew(f, top_global->natoms);

    /* Print to log file  */
    walltime_accounting_start(walltime_accounting);
    wallcycle_start(wcycle, ewcRUN);
    print_start(fplog, cr, walltime_accounting, "Test Particle Insertion");

    /* The last charge group is the group to be inserted */
    cg_tp = top->cgs.nr - 1;
    a_tp0 = top->cgs.index[cg_tp];
    a_tp1 = top->cgs.index[cg_tp+1];
    if (debug)
    {
        fprintf(debug, "TPI cg %d, atoms %d-%d\n", cg_tp, a_tp0, a_tp1);
    }
    if (a_tp1 - a_tp0 > 1 &&
        (inputrec->rlist < inputrec->rcoulomb ||
         inputrec->rlist < inputrec->rvdw))
    {
        gmx_fatal(FARGS, "Can not do TPI for multi-atom molecule with a twin-range cut-off");
    }
    snew(x_mol, a_tp1-a_tp0);

    bDispCorr = (inputrec->eDispCorr != edispcNO);
    bCharge   = FALSE;
    for (i = a_tp0; i < a_tp1; i++)
    {
        /* Copy the coordinates of the molecule to be insterted */
        copy_rvec(state->x[i], x_mol[i-a_tp0]);
        /* Check if we need to print electrostatic energies */
        bCharge |= (mdatoms->chargeA[i] != 0 ||
                    (mdatoms->chargeB && mdatoms->chargeB[i] != 0));
    }
    bRFExcl = (bCharge && EEL_RF(fr->eeltype) && fr->eeltype != eelRF_NEC);

    calc_cgcm(fplog, cg_tp, cg_tp+1, &(top->cgs), state->x, fr->cg_cm);
    if (bCavity)
    {
        if (norm(fr->cg_cm[cg_tp]) > 0.5*inputrec->rlist && fplog)
        {
            fprintf(fplog, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
            fprintf(stderr, "WARNING: Your TPI molecule is not centered at 0,0,0\n");
        }
    }
    else
    {
        /* Center the molecule to be inserted at zero */
        for (i = 0; i < a_tp1-a_tp0; i++)
        {
            rvec_dec(x_mol[i], fr->cg_cm[cg_tp]);
        }
    }

    if (fplog)
    {
        fprintf(fplog, "\nWill insert %d atoms %s partial charges\n",
                a_tp1-a_tp0, bCharge ? "with" : "without");

        fprintf(fplog, "\nWill insert %d times in each frame of %s\n",
                (int)nsteps, opt2fn("-rerun", nfile, fnm));
    }

int write_trxframe_indexed(t_trxstatus *status, const t_trxframe *fr, int nind,
                           const int *ind, gmx_conect gc)
{
    char  title[STRLEN];
    rvec *xout = NULL, *vout = NULL, *fout = NULL;
    int   i, ftp = -1;
    real  prec;

    if (fr->bPrec)
    {
        prec = fr->prec;
    }
    else
    {
        prec = 1000.0;
    }

    if (status->tng)
    {
        ftp = efTNG;
    }
    else if (status->fio)
    {
        ftp = gmx_fio_getftp(status->fio);
    }
    else
    {
        gmx_incons("No input file available");
    }

    switch (ftp)
    {
        case efTRR:
        case efTNG:
            break;
        default:
            if (!fr->bX)
            {
                gmx_fatal(FARGS, "Need coordinates to write a %s trajectory",
                          ftp2ext(ftp));
            }
            break;
    }

    switch (ftp)
    {
        case efTRR:
        case efTNG:
            if (fr->bV)
            {
                snew(vout, nind);
                for (i = 0; i < nind; i++)
                {
                    copy_rvec(fr->v[ind[i]], vout[i]);
                }
            }
            if (fr->bF)
            {
                snew(fout, nind);
                for (i = 0; i < nind; i++)
                {
                    copy_rvec(fr->f[ind[i]], fout[i]);
                }
            }
        /* no break */
        case efXTC:
            if (fr->bX)
            {
                snew(xout, nind);
                for (i = 0; i < nind; i++)
                {
                    copy_rvec(fr->x[ind[i]], xout[i]);
                }
            }
            break;
        default:
            break;
    }

    switch (ftp)
    {
        case efTNG:
            gmx_write_tng_from_trxframe(status->tng, fr, nind);
            break;
        case efXTC:
            write_xtc(status->fio, nind, fr->step, fr->time, fr->box, xout, prec);
            break;
        case efTRR:
            gmx_trr_write_frame(status->fio, nframes_read(status),
                                fr->time, fr->step, fr->box, nind, xout, vout, fout);
            break;
        case efGRO:
        case efPDB:
        case efBRK:
        case efENT:
            if (!fr->bAtoms)
            {
                gmx_fatal(FARGS, "Can not write a %s file without atom names",
                          ftp2ext(ftp));
            }
//.........这里部分代码省略.........

static void mk_diamond(t_atoms *a,rvec x[],real odist,t_symtab *symtab,
		       gmx_bool bPBC,matrix box)
{
  
  int   i,ib,j,k,l,m,nrm=0;
  t_bbb *bbb;
  gmx_bool  *bRemove;
  rvec  dx;
  
  do {
    nrm = 0;
    bbb = mk_bonds(a->nr,x,odist,bPBC,box);
    
    for(i=0; (i<a->nr); i++) {
      if (bbb[i].n < 2) {
	for(k=0; (k<bbb[i].n); k++) {
	  ib = bbb[i].aa[k];
	  for(j=0; (j<bbb[ib].n); j++)
	    if (bbb[ib].aa[j] == i)
	      break;
	  if (j == bbb[ib].n)
	    gmx_fatal(FARGS,"Bond inconsistency (%d not in list of %d)!\n",i,ib);
	  for( ; (j<bbb[ib].n-1); j++)
	    bbb[ib].aa[j] = bbb[ib].aa[j+1];
	  bbb[ib].n--;
	  nrm++;
	}
	bbb[i].n = 0;
      }
    }
  
    for(i=j=0; (i<a->nr); i++) {
      if (bbb[i].n >= 2) {
	copy_rvec(x[i],x[j]);
	j++;
      }
    }
    fprintf(stderr,"Kicking out %d carbon atoms (out of %d)\n",
	    a->nr-j,a->nr);
    a->nr = j;
    sfree(bbb);
  } while (nrm > 0);
  /* Rename atoms */
  bbb = mk_bonds(a->nr,x,odist,bPBC,box);
  for(i=0; (i<a->nr); i++) {
    switch (bbb[i].n) {
    case 4:
      a->atomname[i] = put_symtab(symtab,"C");
      break;
    case 3:
      a->atomname[i] = put_symtab(symtab,"CH1");
      break;
    case 2:
      a->atomname[i] = put_symtab(symtab,"CH2");
      break;
    default:
      gmx_fatal(FARGS,"This atom (%d) has %d bonds only",i,bbb[i].n);
    }
  }
  sfree(bbb);
}

//.........这里部分代码省略.........
		      xvgr_tlabel(),"(degrees)");
    if (bPrintXvgrCodes())
      fprintf(fkinkr,"@ subtitle \"+ = ) (   - = ( )\"\n");
    fkinkl = xvgropen(opt2fn("-okl",NFILE,fnm),"Lateral kink angles",
		      xvgr_tlabel(),"(degrees)");
  }

  if (opt2bSet("-oa",NFILE,fnm)) {
    init_t_atoms(&outatoms,3*n,FALSE);
    outatoms.nr = 3*n;
    for(i=0; i<3*n; i++) {
      outatoms.atomname[i] = &anm;
      outatoms.atom[i].resnr = i/3;
      outatoms.resname[i/3] = &rnm;
    }
    fpdb = open_trx(opt2fn("-oa",NFILE,fnm),"w");
  } else
    fpdb = -1;
  
  read_first_frame(&status,ftp2fn(efTRX,NFILE,fnm),&fr,TRX_NEED_X); 
  
  do {
    rm_pbc(&top.idef,ePBC,fr.natoms,fr.box,fr.x,fr.x);
    calc_axes(fr.x,top.atoms.atom,gnx,index,!bZ,&bun);
    t = convert_time(fr.time);
    fprintf(flen," %10g",t);
    fprintf(fdist," %10g",t);
    fprintf(fz," %10g",t);
    fprintf(ftilt," %10g",t);
    fprintf(ftiltr," %10g",t);
    fprintf(ftiltl," %10g",t);
    if (bKink) {
      fprintf(fkink," %10g",t);
      fprintf(fkinkr," %10g",t);
      fprintf(fkinkl," %10g",t);
    }

    for(i=0; i<bun.n; i++) {
      fprintf(flen," %6g",bun.len[i]);
      fprintf(fdist," %6g",norm(bun.mid[i]));
      fprintf(fz," %6g",bun.mid[i][ZZ]);
      fprintf(ftilt," %6g",RAD2DEG*acos(bun.dir[i][ZZ]));
      comp = bun.mid[i][XX]*bun.dir[i][XX]+bun.mid[i][YY]*bun.dir[i][YY];
      fprintf(ftiltr," %6g",RAD2DEG*
	      asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
      comp = bun.mid[i][YY]*bun.dir[i][XX]-bun.mid[i][XX]*bun.dir[i][YY];
      fprintf(ftiltl," %6g",RAD2DEG*
	      asin(comp/sqrt(sqr(comp)+sqr(bun.dir[i][ZZ]))));
      if (bKink) {
	rvec_sub(bun.end[0][i],bun.end[2][i],va);
	rvec_sub(bun.end[2][i],bun.end[1][i],vb);
	unitv_no_table(va,va);
	unitv_no_table(vb,vb);
	fprintf(fkink," %6g",RAD2DEG*acos(iprod(va,vb)));
	cprod(va,vb,vc);
	copy_rvec(bun.mid[i],vr);
	vr[ZZ] = 0;
	unitv_no_table(vr,vr);
	fprintf(fkinkr," %6g",RAD2DEG*asin(iprod(vc,vr)));
	vl[XX] = vr[YY];
	vl[YY] = -vr[XX];
	vl[ZZ] = 0;
	fprintf(fkinkl," %6g",RAD2DEG*asin(iprod(vc,vl)));
      }
    }
    fprintf(flen,"\n");
    fprintf(fdist,"\n");
    fprintf(fz,"\n");
    fprintf(ftilt,"\n");
    fprintf(ftiltr,"\n");
    fprintf(ftiltl,"\n");
    if (bKink) {
      fprintf(fkink,"\n");
      fprintf(fkinkr,"\n");
      fprintf(fkinkl,"\n");
    }
    if (fpdb >= 0)
      dump_axes(fpdb,&fr,&outatoms,&bun);
  } while(read_next_frame(status,&fr));

  close_trx(status);
  
  if (fpdb >= 0)
    close_trx(fpdb);
  fclose(flen);
  fclose(fdist);
  fclose(fz);
  fclose(ftilt);
  fclose(ftiltr);
  fclose(ftiltl);
  if (bKink) {
    fclose(fkink);
    fclose(fkinkr);
    fclose(fkinkl);
  }
  
  thanx(stderr);
  
  return 0;
}

void read_eigenvectors(const char *file, int *natoms, gmx_bool *bFit,
                       rvec **xref, gmx_bool *bDMR,
                       rvec **xav, gmx_bool *bDMA,
                       int *nvec, int **eignr,
                       rvec ***eigvec, real **eigval)
{
    gmx_trr_header_t   head;
    int                i, snew_size;
    struct t_fileio   *status;
    rvec              *x;
    matrix             box;
    gmx_bool           bOK;

    *bDMR = FALSE;

    /* read (reference (t=-1) and) average (t=0) structure */
    status = gmx_trr_open(file, "r");
    gmx_trr_read_frame_header(status, &head, &bOK);
    *natoms = head.natoms;
    snew(*xav, *natoms);
    gmx_trr_read_frame_data(status, &head, box, *xav, NULL, NULL);

    if ((head.t >= -1.1) && (head.t <= -0.9))
    {
        snew(*xref, *natoms);
        for (i = 0; i < *natoms; i++)
        {
            copy_rvec((*xav)[i], (*xref)[i]);
        }
        *bDMR = (head.lambda > 0.5);
        *bFit = (head.lambda > -0.5);
        if (*bFit)
        {
            fprintf(stderr, "Read %smass weighted reference structure with %d atoms from %s\n", *bDMR ? "" : "non ", *natoms, file);
        }
        else
        {
            fprintf(stderr, "Eigenvectors in %s were determined without fitting\n", file);
            sfree(*xref);
            *xref = NULL;
        }
        gmx_trr_read_frame_header(status, &head, &bOK);
        gmx_trr_read_frame_data(status, &head, box, *xav, NULL, NULL);
    }
    else
    {
        *bFit = TRUE;
        *xref = NULL;
    }
    *bDMA = (head.lambda > 0.5);
    if ((head.t <= -0.01) || (head.t >= 0.01))
    {
        fprintf(stderr, "WARNING: %s does not start with t=0, which should be the "
                "average structure. This might not be a eigenvector file. "
                "Some things might go wrong.\n",
                file);
    }
    else
    {
        fprintf(stderr,
                "Read %smass weighted average/minimum structure with %d atoms from %s\n",
                *bDMA ? "" : "non ", *natoms, file);
    }

    snew(x, *natoms);
    snew_size = 10;
    snew(*eignr, snew_size);
    snew(*eigval, snew_size);
    snew(*eigvec, snew_size);

    *nvec = 0;
    while (gmx_trr_read_frame_header(status, &head, &bOK))
    {
        gmx_trr_read_frame_data(status, &head, box, x, NULL, NULL);
        if (*nvec >= snew_size)
        {
            snew_size += 10;
            srenew(*eignr, snew_size);
            srenew(*eigval, snew_size);
            srenew(*eigvec, snew_size);
        }
        i                = head.step;
        (*eigval)[*nvec] = head.t;
        (*eignr)[*nvec]  = i-1;
        snew((*eigvec)[*nvec], *natoms);
        for (i = 0; i < *natoms; i++)
        {
            copy_rvec(x[i], (*eigvec)[*nvec][i]);
        }
        (*nvec)++;
    }
    sfree(x);
    gmx_trr_close(status);
    fprintf(stderr, "Read %d eigenvectors (for %d atoms)\n\n", *nvec, *natoms);
}

void calc_running_com(t_pull *pull) {
  int i,j,n;
  rvec ave;
  real tm;
  
  /* for group i, we have nhist[i] points of history. The maximum nr of points
     is pull->reflag. The array comhist[i][0..nhist[i]-1] has the positions
     of the center of mass over the last nhist[i] steps. x_unc[i] has the
     new ones. Remove the oldest one from the list, add the new one, calculate
     the average, put that in x_unc instead and return. We just move all coms
     1 down in the list and add the latest one to the top. 
   */
  
  if (pull->bCyl) {
    /* act on dyna groups */
    for (i=0;i<pull->pull.n;i++) {
      clear_rvec(ave);
      for (j=0;j<(pull->reflag-1);j++) {
	copy_rvec(pull->dyna.comhist[i][j+1],pull->dyna.comhist[i][j]);
	rvec_add(ave,pull->dyna.comhist[i][j],ave);
      }
      copy_rvec(pull->dyna.x_unc[i],pull->dyna.comhist[i][j]);
      rvec_add(ave,pull->dyna.x_unc[i],ave);
      svmul(1.0/pull->reflag,ave,ave); 

      /* now ave has the running com for group i, copy it to x_unc[i] */
      copy_rvec(ave,pull->dyna.x_unc[i]);

#ifdef DEBUG
      if (pull->bVerbose) 
	for (n=0;n<pull->reflag;n++) 
	  fprintf(stderr,"Comhist %d, grp %d: %8.3f%8.3f%8.3f\n",n,i,
		  pull->dyna.comhist[i][n][XX],
		  pull->dyna.comhist[i][n][YY],
  		  pull->dyna.comhist[i][n][ZZ]);
#endif
    }
  } else {
    /* act on ref group */
    clear_rvec(ave);
    
    for (j=0;j<(pull->reflag-1);j++) {
      copy_rvec(pull->ref.comhist[0][j+1],pull->ref.comhist[0][j]);
      rvec_add(ave,pull->ref.comhist[0][j],ave);
    }
    
    copy_rvec(pull->ref.x_unc[0],pull->ref.comhist[0][j]);
    rvec_add(ave,pull->ref.x_unc[0],ave);
    svmul(1.0/pull->reflag,ave,ave); 
    /* now ave has the running com for group i, copy it to x_unc[0] */
    copy_rvec(ave,pull->ref.x_unc[0]);

#ifdef DEBUG
    if (pull->bVerbose) 
      for (i=0;i<pull->reflag;i++) 
	fprintf(stderr,"Comhist %d: %8.3f%8.3f%8.3f\n",i,
		pull->ref.comhist[0][i][XX],
		pull->ref.comhist[0][i][YY],
		pull->ref.comhist[0][i][ZZ]);
#endif
  }
}

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];
//.........这里部分代码省略.........

void dd_pmeredist_f(struct gmx_pme_t *pme, pme_atomcomm_t *atc,
                    int n, rvec *f,
                    gmx_bool bAddF)
{
    int *commnode, *buf_index;
    int  nnodes_comm, local_pos, buf_pos, i, scount, rcount, node;

    commnode  = atc->node_dest;
    buf_index = atc->buf_index;

    nnodes_comm = min(2*atc->maxshift, atc->nslab-1);

    local_pos = atc->count[atc->nodeid];
    buf_pos   = 0;
    for (i = 0; i < nnodes_comm; i++)
    {
        scount = atc->rcount[i];
        rcount = atc->count[commnode[i]];
        if (scount > 0 || rcount > 0)
        {
            /* Communicate the forces */
            pme_dd_sendrecv(atc, TRUE, i,
                            atc->f[local_pos], scount*sizeof(rvec),
                            pme->bufv[buf_pos], rcount*sizeof(rvec));
            local_pos += scount;
        }
        buf_index[commnode[i]] = buf_pos;
        buf_pos               += rcount;
    }

    local_pos = 0;
    if (bAddF)
    {
        for (i = 0; i < n; i++)
        {
            node = atc->pd[i];
            if (node == atc->nodeid)
            {
                /* Add from the local force array */
                rvec_inc(f[i], atc->f[local_pos]);
                local_pos++;
            }
            else
            {
                /* Add from the receive buffer */
                rvec_inc(f[i], pme->bufv[buf_index[node]]);
                buf_index[node]++;
            }
        }
    }
    else
    {
        for (i = 0; i < n; i++)
        {
            node = atc->pd[i];
            if (node == atc->nodeid)
            {
                /* Copy from the local force array */
                copy_rvec(atc->f[local_pos], f[i]);
                local_pos++;
            }
            else
            {
                /* Copy from the receive buffer */
                copy_rvec(pme->bufv[buf_index[node]], f[i]);
                buf_index[node]++;
            }
        }
    }
}

//.........这里部分代码省略.........
        fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
                          "Incremental local helix tilt", "Time(ps)", "Tilt (degrees)",
                          oenv);
        fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
                              "Incremental local helix rotation", "Time(ps)",
                              "Rotation (degrees)", oenv);
    }
    else
    {
        fptilt = xvgropen(opt2fn("-otilt", NFILE, fnm),
                          "Cumulative local helix tilt", "Time(ps)", "Tilt (degrees)", oenv);
        fprotation = xvgropen(opt2fn("-orot", NFILE, fnm),
                              "Cumulative local helix rotation", "Time(ps)",
                              "Rotation (degrees)", oenv);
    }

    clear_rvecs(3, unitaxes);
    unitaxes[0][0] = 1;
    unitaxes[1][1] = 1;
    unitaxes[2][2] = 1;

    gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);

    do
    {
        /* initialisation for correct distance calculations */
        set_pbc(&pbc, ePBC, box);
        /* make molecules whole again */
        gmx_rmpbc(gpbc, natoms, box, x);

        /* copy coords to our smaller arrays */
        for (i = 0; i < iCA; i++)
        {
            copy_rvec(x[ind_CA[i]], x_CA[i]);
            if (bSC)
            {
                copy_rvec(x[ind_SC[i]], x_SC[i]);
            }
        }

        for (i = 0; i < iCA-3; i++)
        {
            rvec_sub(x_CA[i+1], x_CA[i], r12[i]);
            rvec_sub(x_CA[i+2], x_CA[i+1], r23[i]);
            rvec_sub(x_CA[i+3], x_CA[i+2], r34[i]);
            rvec_sub(r12[i], r23[i], diff13[i]);
            rvec_sub(r23[i], r34[i], diff24[i]);
            /* calculate helix axis */
            cprod(diff13[i], diff24[i], helixaxis[i]);
            svmul(1.0/norm(helixaxis[i]), helixaxis[i], helixaxis[i]);

            tmp       = cos_angle(diff13[i], diff24[i]);
            twist[i]  = 180.0/M_PI * std::acos( tmp );
            radius[i] = std::sqrt( norm(diff13[i])*norm(diff24[i]) ) / (2.0* (1.0-tmp) );
            rise[i]   = std::abs(iprod(r23[i], helixaxis[i]));

            svmul(radius[i]/norm(diff13[i]), diff13[i], v1);
            svmul(radius[i]/norm(diff24[i]), diff24[i], v2);

            rvec_sub(x_CA[i+1], v1, residueorigin[i+1]);
            rvec_sub(x_CA[i+2], v2, residueorigin[i+2]);
        }
        residueradius[0] = residuetwist[0] = residuerise[0] = 0;

        residueradius[1] = radius[0];
        residuetwist[1]  = twist[0];

//.........这里部分代码省略.........
    if (bCV)
    {
        snew(sumv, fr.natoms);
    }
    if (bCF)
    {
        snew(sumf, fr.natoms);
    }
    nr_xfr = 0;
    nr_vfr = 0;
    nr_ffr = 0;

    if (bCom && bPBC)
    {
        gpbc = gmx_rmpbc_init(&top.idef, ePBC, fr.natoms);
    }

    do
    {
        time = output_env_conv_time(oenv, fr.time);

        if (fr.bX && bNoJump && fr.bBox)
        {
            if (xp)
            {
                remove_jump(fr.box, fr.natoms, xp, fr.x);
            }
            else
            {
                snew(xp, fr.natoms);
            }
            for (i = 0; i < fr.natoms; i++)
            {
                copy_rvec(fr.x[i], xp[i]);
            }
        }

        if (fr.bX && bCom && bPBC)
        {
            gmx_rmpbc_trxfr(gpbc, &fr);
        }

        if (bVD && fr.bV)
        {
            update_histo(isize[0], index[0], fr.v, &nvhisto, &vhisto, binwidth);
        }

        if (bOX && fr.bX)
        {
            print_data(outx, time, fr.x, mass, bCom, ngroups, isize, index, bDim, sffmt);
        }
        if (bOXT && fr.bX)
        {
            frout = fr;
            if (!frout.bAtoms)
            {
                frout.atoms  = &top.atoms;
                frout.bAtoms = TRUE;
            }
            write_trx_x(status_out, &frout, mass, bCom, ngroups, isize, index);
        }
        if (bOV && fr.bV)
        {
            print_data(outv, time, fr.v, mass, bCom, ngroups, isize, index, bDim, sffmt);
        }
        if (bOF && fr.bF)

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

void sgangle_plot_single(const char *fn,const char *afile,const char *dfile, 
			 const char *d1file, const char *d2file,
			 atom_id index1[], int gnx1, char *grpn1,
			 atom_id index2[], int gnx2, char *grpn2,
			 t_topology *top,int ePBC, const output_env_t oenv)
{
  FILE         
    *sg_angle,           /* xvgr file with angles */
    *sg_distance = NULL, /* xvgr file with distances */
    *sg_distance1 = NULL,/* xvgr file with distance between plane and atom */
    *sg_distance2 = NULL;/* xvgr file with distance between plane and atom2 */
  real         
    t,                   /* time */
    angle,               /* cosine of angle between two groups */
    distance,            /* distance between two groups. */
    distance1,           /* distance between plane and one of two atoms */
    distance2;           /* same for second of two atoms */
  t_trxstatus *status;
  int        natoms,teller=0;
  int        i;
  rvec       *x0;   /* coordinates, and coordinates corrected for pb */
  rvec       *xzero;
  matrix     box;        
  char       buf[256];   /* for xvgr title */
  gmx_rmpbc_t  gpbc=NULL;
  

  if ((natoms = read_first_x(oenv,&status,fn,&t,&x0,box)) == 0)
    gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
  
  sprintf(buf,"Angle between %s and %s",grpn1,grpn2);
  sg_angle = xvgropen(afile,buf,"Time (ps)","Cos(angle) ",oenv);
  
  if (dfile) {
    sprintf(buf,"Distance between %s and %s",grpn1,grpn2);
    sg_distance = xvgropen(dfile,buf,"Time (ps)","Distance (nm)",oenv);
  }
  
  if (d1file) {
    sprintf(buf,"Distance between plane and first atom of vector");
    sg_distance1 = xvgropen(d1file,buf,"Time (ps)","Distance (nm)",oenv);
  }
  
  if (d2file) {
    sprintf(buf,"Distance between plane and second atom of vector");
    sg_distance2 = xvgropen(d2file,buf,"Time (ps","Distance (nm)",oenv);
  }
  
  snew(xzero,natoms);
  gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);

  do {
    teller++;
    
    gmx_rmpbc(gpbc,natoms,box,x0);
    if (teller==1) {
      for(i=0;i<natoms;i++)
	copy_rvec(x0[i],xzero[i]);
    }
    
    
    calc_angle_single(ePBC,box,
		      xzero,x0,index1,index2,gnx1,gnx2,&angle,
		      &distance,&distance1,&distance2);
    
    fprintf(sg_angle,"%12g  %12g  %12g\n",t,angle,acos(angle)*180.0/M_PI);
    if (dfile)
      fprintf(sg_distance,"%12g  %12g\n",t,distance);
    if (d1file)
      fprintf(sg_distance1,"%12g  %12g\n",t,distance1);
    if (d2file)
      fprintf(sg_distance2,"%12g  %12g\n",t,distance1);
    
  } while (read_next_x(oenv,status,&t,natoms,x0,box));
  gmx_rmpbc_done(gpbc);
  
  fprintf(stderr,"\n");
  close_trj(status);
  ffclose(sg_angle);
  if (dfile)
    ffclose(sg_distance);
  if (d1file)
    ffclose(sg_distance1);
  if (d2file)
    ffclose(sg_distance2);
  
  sfree(x0);
}

/* calculates center of mass of selection index from all coordinates x */
void pull_calc_coms(t_commrec *cr,
                    struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, double t,
                    rvec x[], rvec *xp)
{
    int          g;
    real         twopi_box = 0;
    pull_comm_t *comm;

    comm = &pull->comm;

    if (comm->rbuf == NULL)
    {
        snew(comm->rbuf, pull->ngroup);
    }
    if (comm->dbuf == NULL)
    {
        snew(comm->dbuf, 3*pull->ngroup);
    }

    if (pull->bRefAt && pull->bSetPBCatoms)
    {
        pull_set_pbcatoms(cr, pull, x, comm->rbuf);

        if (cr != NULL && DOMAINDECOMP(cr))
        {
            /* We can keep these PBC reference coordinates fixed for nstlist
             * steps, since atoms won't jump over PBC.
             * This avoids a global reduction at the next nstlist-1 steps.
             * Note that the exact values of the pbc reference coordinates
             * are irrelevant, as long all atoms in the group are within
             * half a box distance of the reference coordinate.
             */
            pull->bSetPBCatoms = FALSE;
        }
    }

    if (pull->cosdim >= 0)
    {
        int m;

        assert(pull->npbcdim <= DIM);

        for (m = pull->cosdim+1; m < pull->npbcdim; m++)
        {
            if (pbc->box[m][pull->cosdim] != 0)
            {
                gmx_fatal(FARGS, "Can not do cosine weighting for trilinic dimensions");
            }
        }
        twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
    }

    for (g = 0; g < pull->ngroup; g++)
    {
        pull_group_work_t *pgrp;

        pgrp = &pull->group[g];

        if (pgrp->bCalcCOM)
        {
            if (pgrp->epgrppbc != epgrppbcCOS)
            {
                dvec   com, comp;
                double wmass, wwmass;
                rvec   x_pbc = { 0, 0, 0 };
                int    i;

                clear_dvec(com);
                clear_dvec(comp);
                wmass  = 0;
                wwmass = 0;

                if (pgrp->epgrppbc == epgrppbcREFAT)
                {
                    /* Set the pbc atom */
                    copy_rvec(comm->rbuf[g], x_pbc);
                }

                for (i = 0; i < pgrp->nat_loc; i++)
                {
                    int  ii, m;
                    real mass, wm;

                    ii   = pgrp->ind_loc[i];
                    mass = md->massT[ii];
                    if (pgrp->weight_loc == NULL)
                    {
                        wm     = mass;
                        wmass += wm;
                    }
                    else
                    {
                        real w;

                        w       = pgrp->weight_loc[i];
                        wm      = w*mass;
                        wmass  += wm;
                        wwmass += wm*w;
                    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            break;
        case 4: /* two or three tetrahedral hydrogens, e.g. -CH3 */
            for (d = 0; (d < DIM); d++)
            {
                xH1[d] = xAI[d]+distH*sin(alfaH)*sb[d]-distH*cos(alfaH)*sij[d];
                xH2[d] = ( xAI[d]
                           - distH*sin(alfaH)*0.5*sb[d]
                           + distH*sin(alfaH)*s6*sa[d]
                           - distH*cos(alfaH)*sij[d] );
                if (xH3[XX] != NOTSET && xH3[YY] != NOTSET && xH3[ZZ] != NOTSET)
                {
                    xH3[d] = ( xAI[d]
                               - distH*sin(alfaH)*0.5*sb[d]
                               - distH*sin(alfaH)*s6*sa[d]
                               - distH*cos(alfaH)*sij[d] );
                }
            }
            break;
        case 5: /* one tetrahedral hydrogen, e.g. C3CH */
        {
            real center;
            rvec dxc;

            for (d = 0; (d < DIM); d++)
            {
                center = (xAJ[d]+xAK[d]+xAL[d])/3.0;
                dxc[d] = xAI[d]-center;
            }
            center = norm(dxc);
            for (d = 0; (d < DIM); d++)
            {
                xH1[d] = xAI[d]+dxc[d]*distH/center;
            }
            break;
        }
        case 6: /* two tetrahedral hydrogens, e.g. C-CH2-C */
        {
            rvec rBB, rCC1, rCC2, rNN;
            real bb, nn;

            for (d = 0; (d < DIM); d++)
            {
                rBB[d] = xAI[d]-0.5*(xAJ[d]+xAK[d]);
            }
            bb = norm(rBB);

            rvec_sub(xAI, xAJ, rCC1);
            rvec_sub(xAI, xAK, rCC2);
            cprod(rCC1, rCC2, rNN);
            nn = norm(rNN);

            for (d = 0; (d < DIM); d++)
            {
                xH1[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb+
                                       sin(alfaH/2.0)*rNN[d]/nn);
                xH2[d] = xAI[d]+distH*(cos(alfaH/2.0)*rBB[d]/bb-
                                       sin(alfaH/2.0)*rNN[d]/nn);
            }
            break;
        }
        case 7: /* two water hydrogens */
            gen_waterhydrogen(2, xa, xh, l);
            break;
        case 10: /* three water hydrogens */
            gen_waterhydrogen(3, xa, xh, l);
            break;
        case 11: /* four water hydrogens */
            gen_waterhydrogen(4, xa, xh, l);
            break;
        case 8: /* two carboxyl oxygens, -COO- */
            for (d = 0; (d < DIM); d++)
            {
                xH1[d] = xAI[d]-distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
                xH2[d] = xAI[d]+distOM*sin(alfaCOM)*sb[d]-distOM*cos(alfaCOM)*sij[d];
            }
            break;
        case 9:          /* carboxyl oxygens and hydrogen, -COOH */
        {
            rvec xa2[4]; /* i,j,k,l   */

            /* first add two oxygens */
            for (d = 0; (d < DIM); d++)
            {
                xH1[d] = xAI[d]-distO *sin(alfaCO )*sb[d]-distO *cos(alfaCO )*sij[d];
                xH2[d] = xAI[d]+distOA*sin(alfaCOA)*sb[d]-distOA*cos(alfaCOA)*sij[d];
            }

            /* now use rule 2 to add hydrogen to 2nd oxygen */
            copy_rvec(xH2, xa2[0]); /* new i = n' */
            copy_rvec(xAI, xa2[1]); /* new j = i  */
            copy_rvec(xAJ, xa2[2]); /* new k = j  */
            copy_rvec(xAK, xa2[3]); /* new l = k, not used */
            calc_h_pos(2, xa2, (xh+2), l);

            break;
        }
        default:
            gmx_fatal(FARGS, "Invalid argument (%d) for nht in routine genh\n", nht);
    } /* end switch */
}