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

void calc_triclinic_images(matrix box,rvec img[]){    int i;    /* Calculate 3 adjacent images in the xy-plane */    copy_rvec(box[0],img[0]);    copy_rvec(box[1],img[1]);    if (img[1][XX] < 0)        svmul(-1,img[1],img[1]);    rvec_sub(img[1],img[0],img[2]);    /* Get the next 3 in the xy-plane as mirror images */    for(i=0; i<3; i++)        svmul(-1,img[i],img[3+i]);    /* Calculate the first 4 out of xy-plane images */    copy_rvec(box[2],img[6]);    if (img[6][XX] < 0)        svmul(-1,img[6],img[6]);    for(i=0; i<3; i++)        rvec_add(img[6],img[i+1],img[7+i]);    /* Mirror the last 4 from the previous in opposite rotation */    for(i=0; i<4; i++)        svmul(-1,img[6 + (2+i) % 4],img[10+i]);}

static void spread_dum3FAD(rvec xi,rvec xj,rvec xk,                           rvec fi,rvec fj,rvec fk,rvec f,real a,real b){    rvec xij,xjk,xperp,Fpij,Fppp,f1,f2,f3;    real a1,b1,c1,c2,invdij,invdij2,invdp,fproj;    int d;    rvec_sub(xj,xi,xij);    rvec_sub(xk,xj,xjk);    /* 6 flops */    invdij = invsqrt(iprod(xij,xij));    invdij2 = invdij * invdij;    c1 = iprod(xij,xjk) * invdij2;    xperp[XX] = xjk[XX] - c1*xij[XX];    xperp[YY] = xjk[YY] - c1*xij[YY];    xperp[ZZ] = xjk[ZZ] - c1*xij[ZZ];    /* xperp in plane ijk, perp. to ij */    invdp = invsqrt(iprod(xperp,xperp));    a1 = a*invdij;    b1 = b*invdp;    /* 45 flops */    /* a1, b1 and c1 are already calculated in constr_dum3FAD       storing them somewhere will save 45 flops!     */    fproj=iprod(xij  ,f)*invdij2;    svmul(fproj,                     xij,  Fpij); /* proj. f on xij */    svmul(iprod(xperp,f)*invdp*invdp,xperp,Fppp); /* proj. f on xperp */    svmul(b1*fproj,                  xperp,f3);    /* 23 flops */    rvec_sub(f,Fpij,f1);  /* f1 = f - Fpij */    rvec_sub(f1,Fppp,f2); /* f2 = f - Fpij - Fppp */    for (d=0; (d<DIM); d++) {        f1[d]*=a1;        f2[d]*=b1;    }    /* 12 flops */    c2=1+c1;    fi[XX] += f[XX] - f1[XX] + c1*f2[XX] + f3[XX];    fi[YY] += f[YY] - f1[YY] + c1*f2[YY] + f3[YY];    fi[ZZ] += f[ZZ] - f1[ZZ] + c1*f2[ZZ] + f3[ZZ];    fj[XX] +=         f1[XX] - c2*f2[XX] - f3[XX];    fj[YY] +=         f1[YY] - c2*f2[YY] - f3[YY];    fj[ZZ] +=         f1[ZZ] - c2*f2[ZZ] - f3[ZZ];    fk[XX] +=                     f2[XX];    fk[YY] +=                     f2[YY];    fk[ZZ] +=                     f2[ZZ];    /* 30 Flops */    /* TOTAL: 113 flops */}

static void calc_axes(rvec x[],t_atom atom[],int gnx[],atom_id *index[],		      gmx_bool bRot,t_bundle *bun){  int  end,i,div,d;  real *mtot,m;  rvec axis[MAX_ENDS],cent;    snew(mtot,bun->n);  for(end=0; end<bun->nend; end++) {    for(i=0; i<bun->n; i++) {      clear_rvec(bun->end[end][i]);      mtot[i] = 0;    }    div = gnx[end]/bun->n;    for(i=0; i<gnx[end]; i++) {      m = atom[index[end][i]].m;      for(d=0; d<DIM; d++)	bun->end[end][i/div][d] += m*x[index[end][i]][d];      mtot[i/div] += m;    }    clear_rvec(axis[end]);    for(i=0; i<bun->n; i++) {      svmul(1.0/mtot[i],bun->end[end][i],bun->end[end][i]);      rvec_inc(axis[end],bun->end[end][i]);    }    svmul(1.0/bun->n,axis[end],axis[end]);  }  sfree(mtot);  rvec_add(axis[0],axis[1],cent);  svmul(0.5,cent,cent);  /* center the bundle on the origin */  for(end=0; end<bun->nend; end++) {    rvec_dec(axis[end],cent);    for(i=0; i<bun->n; i++)      rvec_dec(bun->end[end][i],cent);  }  if (bRot) {    /* rotate the axis parallel to the z-axis */    rotate_ends(bun,axis[0],YY,ZZ);    rotate_ends(bun,axis[0],XX,ZZ);  }  for(i=0; i<bun->n; i++) {    rvec_add(bun->end[0][i],bun->end[1][i],bun->mid[i]);    svmul(0.5,bun->mid[i],bun->mid[i]);    rvec_sub(bun->end[0][i],bun->end[1][i],bun->dir[i]);    bun->len[i] = norm(bun->dir[i]);    unitv(bun->dir[i],bun->dir[i]);  }}

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

static void average_data(rvec x[],rvec xav[],real *mass,			 int ngrps,int isize[],atom_id **index){  int  g,i,ind,d;  real m,mtot;  rvec tmp;  double sum[DIM];  for(g=0; g<ngrps; g++) {    for(d=0; d<DIM; d++)      sum[d] = 0;    clear_rvec(xav[g]);    mtot = 0;    for(i=0; i<isize[g]; i++) {      ind = index[g][i];      if (mass) {	m = mass[ind];	svmul(m,x[ind],tmp);	for(d=0; d<DIM; d++)	  sum[d] += tmp[d];	mtot += m;      } else	for(d=0; d<DIM; d++)	  sum[d] += x[ind][d];    }    if (mass) {      for(d=0; d<DIM; d++)	xav[g][d] = sum[d]/mtot;    } else {      /* mass=NULL, so these are forces: sum only (do not average) */      for(d=0; d<DIM; d++)	xav[g][d] = sum[d];    }  }}

static void calc_comg(int is,int *coi,int *index,bool bMass,t_atom *atom,		      rvec *x,rvec *x_comg){  int  c,i,d;  rvec xc;  real mtot,m;  if (bMass && atom==NULL)    gmx_fatal(FARGS,"No masses available while mass weighting was requested");  for(c=0; c<is; c++) {    clear_rvec(xc);    mtot = 0;    for(i=coi[c]; i<coi[c+1]; i++) {      if (bMass) {	m = atom[index[i]].m;	for(d=0; d<DIM; d++)	  xc[d] += m*x[index[i]][d];	mtot += m;      } else {	rvec_inc(xc,x[index[i]]);	mtot += 1.0;      }    }    svmul(1/mtot,xc,x_comg[c]);  }}

void unitcell_d(rvec x[],rvec box,real odist){  rvec cc[8] = {    { 0,   0,         0 }, /* C1 */    { cx, cy,       -cz }, /* C2 */    { cx, cy+cy2,     0 }, /* C3 */    { 0,  2*cy+cy2, -cz }, /* C4 */    { 0,   0,         1 }, /* C5 */    { cx, cy,      1+cz }, /* C6 */    { cx, cy+cy2,     1 }, /* C7 */    { 0,  2*cy+cy2,1+cz }, /* C8 */  };  rvec hh[4] = {    { 0,   0,         1  }, /* H relative to C */    { cx,  cy,       -cz },    { cx, -cy,       -cz },     {-cy2,  0,       -cz }  };  int  i,iin,iout,j,m;  rvec tmp,t2,dip;    clear_rvec(dip);  for(i=0; (i<8); i++) {    iin  = i;    iout = i;    svmul(odist,cc[iin],x[iout]);  }    box[XX] = 2*cx;  box[YY] = 2*(cy2+cy);  box[ZZ] = 2*(1+cz);  for(i=0; (i<DIM); i++)    box[i] *= odist;      printf("Unitcell:  %10.5f  %10.5f  %10.5f/n",box[XX],box[YY],box[ZZ]);}

static void decrease_step_size(int nshell,t_shell s[]){  int i;    for(i=0; i<nshell; i++)    svmul(0.8,s[i].step,s[i].step);}

/* Determine the (weighted) sum vector from positions x */extern double get_sum_of_positions(rvec x[], real weight[], const int nat, dvec dsumvec){    int i;    rvec x_weighted;    double weight_sum = 0.0;    /* Zero out the center */    clear_dvec(dsumvec);    /* Loop over all atoms and add their weighted position vectors */    if (weight != NULL)    {        for (i=0; i<nat; i++)        {            weight_sum += weight[i];            svmul(weight[i], x[i], x_weighted);            dsumvec[XX] += x_weighted[XX];            dsumvec[YY] += x_weighted[YY];            dsumvec[ZZ] += x_weighted[ZZ];        }    }    else    {        for (i=0; i<nat; i++)        {            dsumvec[XX] += x[i][XX];            dsumvec[YY] += x[i][YY];            dsumvec[ZZ] += x[i][ZZ];        }    }    return weight_sum;}

/*! * /param[in]  top   Topology structure with masses. * /param[in]  f     Forces on all atoms. * /param[in]  block t_block structure that divides /p index into blocks. * /param[in]  index Indices of atoms. * /param[out] fout  /p block->nr Forces on COG positions. */voidgmx_calc_cog_f_block(const gmx_mtop_t *top, rvec f[], const t_block *block, const int index[],                     rvec fout[]){    GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),                       "No masses available while mass weighting was requested");    int molb = 0;    for (int b = 0; b < block->nr; ++b)    {        rvec fb;        clear_rvec(fb);        real mtot = 0;        for (int i = block->index[b]; i < block->index[b+1]; ++i)        {            const int  ai   = index[i];            const real mass = mtopGetAtomMass(top, ai, &molb);            for (int d = 0; d < DIM; ++d)            {                fb[d] += f[ai][d] / mass;            }            mtot += mass;        }        svmul(mtot / (block->index[b+1] - block->index[b]), fb, fout[b]);    }}

void make_refgrps(t_pull *pull,matrix box,t_mdatoms *md) {  int ngrps,i,ii,j,k,m;  static bool bFirst = TRUE;  real dr,mass;  real truemass;  rvec test;  ngrps = pull->pull.n;  if (bFirst) {    snew(pull->dyna.ngx,ngrps);    snew(pull->dyna.idx,ngrps);    snew(pull->dyna.weights,ngrps);    for (i=0;i<ngrps;i++) {      snew(pull->dyna.idx[i],pull->ref.ngx[0]);    /* more than nessary */      snew(pull->dyna.weights[i],pull->ref.ngx[0]);    }          bFirst = FALSE;  }  /* loop over all groups to make a reference group for each*/  for (i=0;i<ngrps;i++) {    k=0;    truemass=0;    pull->dyna.tmass[i] = 0;    pull->dyna.ngx[i] = 0;    /* loop over all atoms in the main ref group */    for (j=0;j<pull->ref.ngx[0];j++) {      ii = pull->ref.idx[0][j];      /* get_distance takes pbc into account */      dr = get_cylinder_distance(pull->ref.x0[0][j],pull->pull.x_unc[i],box);      if (dr < pull->rc) {	/* add to index, to sum of COM, to weight array */	mass = md->massT[ii];	truemass += mass;	pull->dyna.ngx[i]++;	pull->dyna.weights[i][k] = get_weight(dr,pull->r,pull->rc);	pull->dyna.idx[i][k] = ii;	for (m=0;m<DIM;m++) 	  pull->dyna.x_unc[i][m] += mass*pull->dyna.weights[i][k]*	    pull->ref.x0[0][j][m];	pull->dyna.tmass[i] += mass*pull->dyna.weights[i][k];	k++;      }    }    /* normalize the new 'x_unc' */    svmul(1/pull->dyna.tmass[i],pull->dyna.x_unc[i],pull->dyna.x_unc[i]);    if (pull->bVerbose)       fprintf(stderr,"Made group %d:%8.3f%8.3f%8.3f wm:%8.3f m:%8.3f/n",	      i,pull->dyna.x_unc[i][0],pull->dyna.x_unc[i][1],	      pull->dyna.x_unc[i][2],pull->dyna.tmass[i],truemass);  }}

real calc_geom(int isize, atom_id *index, rvec *x, rvec geom_center, rvec min,               rvec max, gmx_bool bDiam){    real diam2, d;    char *grpnames;    int ii, i, j;    clear_rvec(geom_center);    diam2 = 0;    if (isize == 0)    {        clear_rvec(min);        clear_rvec(max);    }    else    {        if (index)            ii = index[0];        else            ii = 0;        for (j = 0; j < DIM; j++)            min[j] = max[j] = x[ii][j];        for (i = 0; i < isize; i++)        {            if (index)                ii = index[i];            else                ii = i;            rvec_inc(geom_center, x[ii]);            for (j = 0; j < DIM; j++)            {                if (x[ii][j] < min[j])                    min[j] = x[ii][j];                if (x[ii][j] > max[j])                    max[j] = x[ii][j];            }            if (bDiam)            {                if (index)                    for (j = i + 1; j < isize; j++)                    {                        d = distance2(x[ii], x[index[j]]);                        diam2 = max(d,diam2);                    }                else                    for (j = i + 1; j < isize; j++)                    {                        d = distance2(x[i], x[j]);                        diam2 = max(d,diam2);                    }            }        }        svmul(1.0 / isize, geom_center, geom_center);    }    return sqrt(diam2);}

void put_atoms_in_triclinic_unitcell(int ecenter, matrix box,                                     int natoms, rvec x[]){    rvec   box_center, shift_center;    real   shm01, shm02, shm12, shift;    int    i, m, d;    calc_box_center(ecenter, box, box_center);    /* The product of matrix shm with a coordinate gives the shift vector       which is required determine the periodic cell position */    shm01 = box[1][0]/box[1][1];    shm02 = (box[1][1]*box[2][0] - box[2][1]*box[1][0])/(box[1][1]*box[2][2]);    shm12 = box[2][1]/box[2][2];    clear_rvec(shift_center);    for (d = 0; d < DIM; d++)    {        rvec_inc(shift_center, box[d]);    }    svmul(0.5, shift_center, shift_center);    rvec_sub(box_center, shift_center, shift_center);    shift_center[0] = shm01*shift_center[1] + shm02*shift_center[2];    shift_center[1] = shm12*shift_center[2];    shift_center[2] = 0;    for (i = 0; (i < natoms); i++)    {        for (m = DIM-1; m >= 0; m--)        {            shift = shift_center[m];            if (m == 0)            {                shift += shm01*x[i][1] + shm02*x[i][2];            }            else if (m == 1)            {                shift += shm12*x[i][2];            }            while (x[i][m]-shift < 0)            {                for (d = 0; d <= m; d++)                {                    x[i][d] += box[m][d];                }            }            while (x[i][m]-shift >= box[m][m])            {                for (d = 0; d <= m; d++)                {                    x[i][d] -= box[m][d];                }            }        }    }}

static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,                         atom_id index[], rvec xref, int ePBC){    const real tol = 1e-4;    gmx_bool   bChanged;    int        m, j, ai, iter;    real       mass, mtot;    rvec       dx, xtest;    /* First simple calculation */    clear_rvec(xref);    mtot = 0;    for (m = 0; (m < nrefat); m++)    {        ai   = index[m];        mass = top->atoms.atom[ai].m;        for (j = 0; (j < DIM); j++)        {            xref[j] += mass*x[ai][j];        }        mtot += mass;    }    svmul(1/mtot, xref, xref);    /* Now check if any atom is more than half the box from the COM */    if (ePBC != epbcNONE)    {        iter = 0;        do        {            bChanged = FALSE;            for (m = 0; (m < nrefat); m++)            {                ai   = index[m];                mass = top->atoms.atom[ai].m/mtot;                pbc_dx(pbc, x[ai], xref, dx);                rvec_add(xref, dx, xtest);                for (j = 0; (j < DIM); j++)                {                    if (std::abs(xtest[j]-x[ai][j]) > tol)                    {                        /* Here we have used the wrong image for contributing to the COM */                        xref[j] += mass*(xtest[j]-x[ai][j]);                        x[ai][j] = xtest[j];                        bChanged = TRUE;                    }                }            }            if (bChanged)            {                printf("COM: %8.3f  %8.3f  %8.3f  iter = %d/n", xref[XX], xref[YY], xref[ZZ], iter);            }            iter++;        }        while (bChanged);    }}

PyObject *apply_rotation( PyObject *self, PyObject *args){  PyObject *Rotation;  PyObject *py_v;  real phi;  if(!PyArg_ParseTuple(args,"OOd",&Rotation,&py_v, &phi))    return NULL;  PyObject *rm1 = PyObject_GetAttrString(Rotation,"m1");  PyObject *rm2 = PyObject_GetAttrString(Rotation,"m2");  matrix m1, m2;  PyObject2matrix(rm1, m1);  PyObject2matrix(rm2, m2);  rvec v, v2;  PyObject *py_v2 = PyObject_GetAttrString(Rotation,"v2");  Pyvec2rvec(py_v, v);  Pyvec2rvec(py_v2, v2);  rvec vec;  rvec_sub(v, v2, vec );  rvec b, d, a, c, e;  mvmul(m1, vec, b);  mvmul(m2, vec, d);  real cc = cos(phi);  svmul( cc, vec, a);  svmul( -cc, b, c);  svmul( sin(phi), d, e);  clear_rvec( vec );  rvec_add( a, b, vec);  rvec_add( vec, c, vec );  rvec_add( vec, e, vec );    clear_rvec(v);  rvec_add( v2, vec, v);  return Py_BuildValue("[ddd]", v[XX], v[YY], v[ZZ] ); }

static void scale_velocities(t_state *state, real fac){    int i;    if (as_rvec_array(state->v.data()))    {        for (i = 0; i < state->natoms; i++)        {            svmul(fac, state->v[i], state->v[i]);        }    }}

static void scale_velocities(t_state *state, real fac){    int i;    if (state->v)    {        for (i = 0; i < state->natoms; i++)        {            svmul(fac, state->v[i], state->v[i]);        }    }}

//! Helper method to calculate a vector from two or three positions.static voidcalc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout){    switch (natoms)    {        case 2:            if (pbc)            {                pbc_dx(pbc, x[1], x[0], xout);            }            else            {                rvec_sub(x[1], x[0], xout);            }            svmul(0.5, xout, cout);            rvec_add(x[0], cout, cout);            break;        case 3:        {            rvec v1, v2;            if (pbc)            {                pbc_dx(pbc, x[1], x[0], v1);                pbc_dx(pbc, x[2], x[0], v2);            }            else            {                rvec_sub(x[1], x[0], v1);                rvec_sub(x[2], x[0], v2);            }            cprod(v1, v2, xout);            rvec_add(x[0], x[1], cout);            rvec_add(cout, x[2], cout);            svmul(1.0/3.0, cout, cout);            break;        }        default:            GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");    }}

/*! * /param[in]  top    Topology structure with masses. * /param[in]  x      Position vectors of all atoms. * /param[in]  pbc    Periodic boundary conditions structure. * /param[in]  nrefat Number of atoms in the index. * /param[in]  index  Indices of atoms. * /param[out] xout   COM position for the indexed atoms. * * Works as gmx_calc_com(), but takes into account periodic boundary * conditions: If any atom is more than half the box from the COM, * it is wrapped around and a new COM is calculated. This is repeated * until no atoms violate the condition. * * Modified from src/tools/gmx_sorient.c in Gromacs distribution. */voidgmx_calc_com_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,                 int nrefat, const int index[], rvec xout){    GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),                       "No masses available while mass weighting was requested");    /* First simple calculation */    clear_rvec(xout);    real mtot = 0;    int  molb = 0;    for (int m = 0; m < nrefat; ++m)    {        const int  ai   = index[m];        const real mass = mtopGetAtomMass(top, ai, &molb);        for (int j = 0; j < DIM; ++j)        {            xout[j] += mass * x[ai][j];        }        mtot += mass;    }    svmul(1.0/mtot, xout, xout);    /* Now check if any atom is more than half the box from the COM */    if (pbc)    {        const real tol  = 1e-4;        bool       bChanged;        do        {            bChanged = false;            molb     = 0;            for (int m = 0; m < nrefat; ++m)            {                rvec       dx, xtest;                const int  ai   = index[m];                const real mass = mtopGetAtomMass(top, ai, &molb) / mtot;                pbc_dx(pbc, x[ai], xout, dx);                rvec_add(xout, dx, xtest);                for (int j = 0; j < DIM; ++j)                {                    if (fabs(xtest[j] - x[ai][j]) > tol)                    {                        /* Here we have used the wrong image for contributing to the COM */                        xout[j] += mass * (xtest[j] - x[ai][j]);                        x[ai][j] = xtest[j];                        bChanged = true;                    }                }            }        }        while (bChanged);    }}

void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees){    rvec angle;    svmul(DEG2RAD, angleInDegrees, angle);    box[XX][XX] = vec[XX];    box[YY][XX] = vec[YY]*cos(angle[ZZ]);    box[YY][YY] = vec[YY]*sin(angle[ZZ]);    box[ZZ][XX] = vec[ZZ]*cos(angle[YY]);    box[ZZ][YY] = vec[ZZ]        *(cos(angle[XX])-cos(angle[YY])*cos(angle[ZZ]))/sin(angle[ZZ]);    box[ZZ][ZZ] = sqrt(gmx::square(vec[ZZ])                       -box[ZZ][XX]*box[ZZ][XX]-box[ZZ][YY]*box[ZZ][YY]);}

static void calc_ringh(rvec xattach, rvec xb, rvec xc, rvec xh){    rvec tab, tac;    real n;    /* Add a proton on a ring to atom attach at distance 0.1 nm */    rvec_sub(xattach, xb, tab);    rvec_sub(xattach, xc, tac);    rvec_add(tab, tac, xh);    n = 0.1/norm(xh);    svmul(n, xh, xh);    rvec_inc(xh, xattach);}

voidgmx_calc_cog(const gmx_mtop_t * /* top */, rvec x[], int nrefat, const int index[], rvec xout){    int                 m, ai;    clear_rvec(xout);    for (m = 0; m < nrefat; ++m)    {        ai = index[m];        rvec_inc(xout, x[ai]);    }    svmul(1.0/nrefat, xout, xout);}

void rm_res_pbc(t_atoms *atoms, rvec *x, matrix box){    int  i, start, n, d, nat;    rvec xcg;    start = 0;    nat   = 0;    clear_rvec(xcg);    for (n = 0; n < atoms->nr; n++)    {        if (!is_hydrogen(*(atoms->atomname[n])))        {            nat++;            rvec_inc(xcg, x[n]);        }        if ( (n+1 == atoms->nr) ||             (atoms->atom[n+1].resind != atoms->atom[n].resind) )        {            /* if nat==0 we have only hydrogens in the solvent,               we take last coordinate as cg */            if (nat == 0)            {                nat = 1;                copy_rvec(x[n], xcg);            }            svmul(1.0/nat, xcg, xcg);            for (d = 0; d < DIM; d++)            {                while (xcg[d] < 0)                {                    for (i = start; i <= n; i++)                    {                        x[i][d] += box[d][d];                    }                    xcg[d] += box[d][d];                }                while (xcg[d] >= box[d][d])                {                    for (i = start; i <= n; i++)                    {                        x[i][d] -= box[d][d];                    }                    xcg[d] -= box[d][d];                }            }            start = n+1;            nat   = 0;            clear_rvec(xcg);        }    }}

static void center_molecule(int atomCount, rvec x[]){    rvec center;    clear_rvec(center);    for (int i = 0; i < atomCount; ++i)    {        rvec_inc(center, x[i]);    }    svmul(1.0/atomCount, center, center);    for (int i = 0; i < atomCount; ++i)    {        rvec_dec(x[i], center);    }}

static voidcalc_vec(int natoms, rvec x[], t_pbc *pbc, rvec xout, rvec cout){    switch (natoms)    {        case 2:            if (pbc)            {                pbc_dx(pbc, x[1], x[0], xout);            }            else            {                rvec_sub(x[1], x[0], xout);            }            svmul(0.5, xout, cout);            rvec_add(x[0], cout, cout);            break;        case 3: {            rvec v1, v2;            if (pbc)            {                pbc_dx(pbc, x[1], x[0], v1);                pbc_dx(pbc, x[2], x[0], v2);            }            else            {                rvec_sub(x[1], x[0], v1);                rvec_sub(x[2], x[0], v2);            }            cprod(v1, v2, xout);            rvec_add(x[0], x[1], cout);            rvec_add(cout, x[2], cout);            svmul(1.0/3.0, cout, cout);            break;        }    }}

voidgmx_calc_cog_block(const gmx_mtop_t * /* top */, rvec x[], const t_block *block, const int index[],                   rvec xout[]){    int                 b, i, ai;    rvec                xb;    for (b = 0; b < block->nr; ++b)    {        clear_rvec(xb);        for (i = block->index[b]; i < block->index[b+1]; ++i)        {            ai = index[i];            rvec_inc(xb, x[ai]);        }        svmul(1.0/(block->index[b+1] - block->index[b]), xb, xout[b]);    }}

void init_resize(t_block *ins_at,rvec *r_ins,pos_ins_t *pos_ins,mem_t *mem_p,rvec *r, gmx_bool bALLOW_ASYMMETRY){	int i,j,at,c,outsidesum,gctr=0;    int idxsum=0;    /*sanity check*/    for (i=0;i<pos_ins->pieces;i++)          idxsum+=pos_ins->nidx[i];    if (idxsum!=ins_at->nr)          gmx_fatal(FARGS,"Piecewise sum of inserted atoms not same as size of group selected to insert.");    snew(pos_ins->geom_cent,pos_ins->pieces);    for (i=0;i<pos_ins->pieces;i++)    {    	c=0;    	outsidesum=0;    	for(j=0;j<DIM;j++)    		pos_ins->geom_cent[i][j]=0;    	for(j=0;j<DIM;j++)    		pos_ins->geom_cent[i][j]=0;    	for (j=0;j<pos_ins->nidx[i];j++)    	{    		at=pos_ins->subindex[i][j];    		copy_rvec(r[at],r_ins[gctr]);    		if( (r_ins[gctr][ZZ]<mem_p->zmax) && (r_ins[gctr][ZZ]>mem_p->zmin) )    		{    			rvec_inc(pos_ins->geom_cent[i],r_ins[gctr]);    			c++;    		}    		else    			outsidesum++;    		gctr++;    	}    	if (c>0)    		svmul(1/(double)c,pos_ins->geom_cent[i],pos_ins->geom_cent[i]);    	if (!bALLOW_ASYMMETRY)    		pos_ins->geom_cent[i][ZZ]=mem_p->zmed;    	fprintf(stderr,"Embedding piece %d with center of geometry: %f %f %f/n",i,pos_ins->geom_cent[i][XX],pos_ins->geom_cent[i][YY],pos_ins->geom_cent[i][ZZ]);    }    fprintf(stderr,"/n");}

/*! * /param[in]  top    Topology structure with masses. * /param[in]  x      Position vectors of all atoms. * /param[in]  nrefat Number of atoms in the index. * /param[in]  index  Indices of atoms. * /param[out] xout   COM position for the indexed atoms. * * Works exactly as gmx_calc_cog() with the exception that a center of * mass are calculated, and hence a topology with masses is required. */voidgmx_calc_com(const gmx_mtop_t *top, rvec x[], int nrefat, const int index[], rvec xout){    GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),                       "No masses available while mass weighting was requested");    clear_rvec(xout);    real mtot = 0;    int  molb = 0;    for (int m = 0; m < nrefat; ++m)    {        const int  ai   = index[m];        const real mass = mtopGetAtomMass(top, ai, &molb);        for (int j = 0; j < DIM; ++j)        {            xout[j] += mass * x[ai][j];        }        mtot += mass;    }    svmul(1.0/mtot, xout, xout);}

/* calculate the com of molecules in x and put it into xa */static void calc_mol_com(int nmol,int *molindex,t_block *mols,t_atoms *atoms,			 rvec *x,rvec *xa){  int  m,mol,i,d;  rvec xm;  real mass,mtot;  for(m=0; m<nmol; m++) {    mol = molindex[m];    clear_rvec(xm);    mtot = 0;    for(i=mols->index[mol]; i<mols->index[mol+1]; i++) {      mass = atoms->atom[i].m;      for(d=0; d<DIM; d++)	xm[d] += mass*x[i][d];      mtot += mass;    }    svmul(1/mtot,xm,xa[m]);  }}