C++ rvec_add函数代码示例

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

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

PyObject *wrap_fit(PyObject *self,PyObject *args)
{

  PyObject *cs1, *cs2, *mass;
  if(!PyArg_ParseTuple(args,"OOO",&cs1, &cs2, &mass))
    return NULL;
  int natoms1 = PySequence_Length(cs1);
  int natoms2 = PySequence_Length(cs2);
  if( natoms1 != natoms2 ) {
    Error("Cannot fit coordinate sets with different lengths");
  }
  rvec x1[natoms1];
  rvec x2[natoms1];
  real m[natoms1];
  PyObject2rvec( cs1, x1, natoms1);
  PyObject2rvec( cs2, x2, natoms2);
  PyObject2real_array(mass, m, natoms1);
  rvec cent;
  center_and_get_vec(x1, natoms1, cent);     // center x1 and get vector for back translation
  center(x2, natoms1);                                // center x2
  do_fit(natoms1, m, x1, x2);               

  int i;
  for(i=0;i<natoms1;i++)                    // translate back
    rvec_add( x2[i], cent, x2[i]);

  PyObject *ret = rvec2PyObject(x2, natoms1);
  return ret;
}

static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
                          real *rmin, real *rmax, int *min_ind)
{
#define NSHIFT 26
    int  sx, sy, sz, i, j, s;
    real sqr_box, r2min, r2max, r2;
    rvec shift[NSHIFT], d0, d;

    sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ]))));

    s = 0;
    for (sz = -1; sz <= 1; sz++)
    {
        for (sy = -1; sy <= 1; sy++)
        {
            for (sx = -1; sx <= 1; sx++)
            {
                if (sx != 0 || sy != 0 || sz != 0)
                {
                    for (i = 0; i < DIM; i++)
                    {
                        shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
                    }
                    s++;
                }
            }
        }
    }

    r2min = sqr_box;
    r2max = 0;

    for (i = 0; i < n; i++)
    {
        for (j = i+1; j < n; j++)
        {
            rvec_sub(x[index[i]], x[index[j]], d0);
            r2 = norm2(d0);
            if (r2 > r2max)
            {
                r2max = r2;
            }
            for (s = 0; s < NSHIFT; s++)
            {
                rvec_add(d0, shift[s], d);
                r2 = norm2(d);
                if (r2 < r2min)
                {
                    r2min      = r2;
                    min_ind[0] = i;
                    min_ind[1] = j;
                }
            }
        }
    }

    *rmin = sqrt(r2min);
    *rmax = sqrt(r2max);
}

void visualize_box(FILE *out,int a0,int r0,matrix box,rvec gridsize)
{
  int     *edge;
  rvec    *vert,shift;
  int     nx,ny,nz,nbox,nat;
  int     i,j,x,y,z;
  int     rectedge[24] = { 0,1, 1,3, 3,2, 0,2, 0,4, 1,5, 3,7, 2,6, 4,5, 5,7, 7,6, 6,4 };

  a0++;
  r0++;
  
  nx = (int)(gridsize[XX]+0.5);
  ny = (int)(gridsize[YY]+0.5);
  nz = (int)(gridsize[ZZ]+0.5);
  nbox = nx*ny*nz;
  if (TRICLINIC(box)) {
    nat = nbox*NCUCVERT;
    snew(vert,nat);
    calc_compact_unitcell_vertices(ecenterDEF,box,vert);
    j = 0;
    for(z=0; z<nz; z++)
      for(y=0; y<ny; y++)
	for(x=0; x<nx; x++) {
	  for(i=0; i<DIM; i++)
	    shift[i] = x*box[0][i]+y*box[1][i]+z*box[2][i];
	  for(i=0; i<NCUCVERT; i++) {
	    rvec_add(vert[i],shift,vert[j]);
	    j++;
	  }
	}
    
    for(i=0; i<nat; i++) {
      fprintf(out,pdbformat,"ATOM",a0+i,"C","BOX",'K'+i/NCUCVERT,r0+i,
	      10*vert[i][XX],10*vert[i][YY],10*vert[i][ZZ]);
      fprintf(out,"\n");
    }
    
    edge = compact_unitcell_edges();
    for(j=0; j<nbox; j++)
      for(i=0; i<NCUCEDGE; i++)
	fprintf(out,"CONECT%5d%5d\n",
		a0 + j*NCUCVERT + edge[2*i],
		a0 + j*NCUCVERT + edge[2*i+1]);
    
    sfree(vert);
  } else {
    i=0;
    for(z=0; z<=1; z++)
      for(y=0; y<=1; y++)
	for(x=0; x<=1; x++) {
	  fprintf(out,pdbformat,"ATOM",a0+i,"C","BOX",'K'+i/8,r0+i,
		  x*10*box[XX][XX],y*10*box[YY][YY],z*10*box[ZZ][ZZ]);
	  fprintf(out,"\n");
	  i++;
	}
    for(i=0; i<24; i+=2)
      fprintf(out,"CONECT%5d%5d\n",a0+rectedge[i],a0+rectedge[i+1]);
  }
}

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

void write_conf_to_grofile(const System& system,
                           const std::string& path,
                           const real sigma,
                           const OutputMode mode)
{
    constexpr char ATOM_NAME[2] = "C";
    constexpr char RESIDUE_NAME[4] = "SOL";

    auto open_mode = std::ios::out;

    if (mode == OutputMode::Append)
    {
        open_mode = open_mode | std::ios::app;
    }

    std::ofstream out { path, open_mode };

    out << system.title << '\n'
        << system.num_atoms() << '\n';

    out.setf(std::ios::fixed);
    out.precision(3);

    uint64_t n = 1;

    for (const auto &list : system.cell_lists)
    {
        for (unsigned i = 0; i < list.num_atoms(); ++i)
        {
            const auto x0 = RVec {
                list.xs.at(i * NDIM + XX),
                list.xs.at(i * NDIM + YY),
                list.xs.at(i * NDIM + ZZ)
            };
            const auto xabs = rvec_add(x0, list.origin);

            out << std::setw(5) << std::right << n
                << std::setw(5) << std::left << RESIDUE_NAME
                << std::setw(5) << std::right << ATOM_NAME
                << std::setw(5) << n
                << std::setw(8) << xabs[XX] * sigma
                << std::setw(8) << xabs[YY] * sigma
                << std::setw(8) << xabs[ZZ] * sigma
                << '\n';

            ++n;
        }
    }

    out << std::setw(9) << std::right << system.box_size[0] * sigma << ' '
        << std::setw(9) << system.box_size[1] * sigma << ' '
        << std::setw(9) << system.box_size[2] * sigma << '\n';
}

//! Helper method to calculate a vector from two or three positions.
static void
calc_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.
 */
void
gmx_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);
    }
}

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

static void
calc_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;
        }
    }
}

void dump_pbc(FILE *fp, t_pbc *pbc)
{
    rvec sum_box;

    fprintf(fp, "ePBCDX = %d\n", pbc->ePBCDX);
    pr_rvecs(fp, 0, "box", pbc->box, DIM);
    pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
    pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
    pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
    rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
    pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
    fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
    fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
    if (pbc->ntric_vec > 0)
    {
        pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
        pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
    }
}

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

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