C++ pbc_dx函数代码示例

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

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
Distance::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                       TrajectoryAnalysisModuleData *pdata)
{
    AnalysisDataHandle       dh   = pdata->dataHandle(data_);
    const Selection         &sel1 = pdata->parallelSelection(sel_[0]);
    const Selection         &sel2 = pdata->parallelSelection(sel_[1]);
    rvec                     dx;
    real                     r;
    const SelectionPosition &p1 = sel1.position(0);
    const SelectionPosition &p2 = sel2.position(0);

    if (pbc != NULL)
    {
        pbc_dx(pbc, p1.x(), p2.x(), dx);
    }
    else
    {
        rvec_sub(p1.x(), p2.x(), dx);
    }
    r = norm(dx);
    dh.startFrame(frnr, fr.time);
    dh.setPoint(0, r);
    dh.setPoints(1, 3, dx);
    dh.finishFrame();
}

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

static real dist2(t_pbc *pbc,rvec x,rvec y)
{
  rvec dx;
  
  pbc_dx(pbc,x,y,dx);
  
  return norm2(dx);
}

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

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

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

// **************** PERIODIC BOUNDARY CONDITIONS ****************
//  Get the PBC-aware distance vector between two positions
cvm::rvector colvarproxy_gromacs::position_distance (cvm::atom_pos const &pos1,
				cvm::atom_pos const &pos2) {
  rvec r1, r2, dr;
  r1[0] = pos1.x;
  r1[1] = pos1.y;
  r1[2] = pos1.z;
  r2[0] = pos2.x;
  r2[1] = pos2.y;
  r2[2] = pos2.z;

  pbc_dx(&gmx_pbc, r1, r2, dr);
  return cvm::atom_pos( dr[0], dr[1], dr[2] );
}

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

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 real calc_dist(t_pbc *pbc,rvec x[],t_block *cgs,int icg,int jcg)
{
  int  i,j;
  rvec dx;
  real d2,mindist2=1000;

  for(i=cgs->index[icg]; (i<cgs->index[icg+1]); i++) {
    for(j=cgs->index[jcg]; (j<cgs->index[jcg+1]); j++) {
      pbc_dx(pbc,x[i],x[j],dx);
      d2 = norm2(dx);
      if (d2 < mindist2)
	mindist2 = d2;
    }
  }
  return sqrt(mindist2);
}

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