C++ set_pbc函数代码示例

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

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

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

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

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

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

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

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

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

        pcrd  = &pull->coord[c];

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

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

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

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

    finish_pull(pull_work);
}

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

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

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

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

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

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

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

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

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

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


    calc_box_center(ecenterRECT, box, center);

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

    clear_rvec(tmp);


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

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

        rvec_inc(mj, mt1);

    }

}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  // Override periodic boundary conditions

  RANK_TO_INDEX( world_rank, px,py,pz );

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

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

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

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

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

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

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

    set_pbc(&pbc, ePBC, box);

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

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

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

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

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

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