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

static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp,                             t_mdatoms *md, rvec *x,                             rvec x_pbc){    int a, m;    if (cr && PAR(cr))    {        if (DOMAINDECOMP(cr))        {            if (!ga2la_get_home(cr->dd->ga2la, pgrp->pbcatom, &a))            {                a = -1;            }        }        else        {            a = pgrp->pbcatom;        }        if (a >= 0 && a < md->homenr)        {            copy_rvec(x[a], x_pbc);        }        else        {            clear_rvec(x_pbc);        }    }    else    {        copy_rvec(x[pgrp->pbcatom], x_pbc);    }}

void make_local_shells(t_commrec *cr,t_mdatoms *md,		       struct gmx_shellfc *shfc){  t_shell *shell;  int a0,a1,*ind,nshell,i;  gmx_domdec_t *dd=NULL;  if (PAR(cr)) {    if (DOMAINDECOMP(cr)) {      dd = cr->dd;      a0 = 0;      a1 = dd->nat_home;    } else {      pd_at_range(cr,&a0,&a1);    }  } else {    /* Single node: we need all shells, just copy the pointer */    shfc->nshell = shfc->nshell_gl;    shfc->shell  = shfc->shell_gl;        return;  }  ind = shfc->shell_index_gl;  nshell = 0;  shell  = shfc->shell;   for(i=a0; i<a1; i++) {    if (md->ptype[i] == eptShell) {      if (nshell+1 > shfc->shell_nalloc) {	shfc->shell_nalloc = over_alloc_dd(nshell+1);	srenew(shell,shfc->shell_nalloc);      }      if (dd) {	shell[nshell] = shfc->shell_gl[ind[dd->gatindex[i]]];      } else {	shell[nshell] = shfc->shell_gl[ind[i]];      }      /* With inter-cg shells we can no do shell prediction,       * so we do not need the nuclei numbers.       */      if (!shfc->bInterCG) {	shell[nshell].nucl1   = i + shell[nshell].nucl1 - shell[nshell].shell;	if (shell[nshell].nnucl > 1)	  shell[nshell].nucl2 = i + shell[nshell].nucl2 - shell[nshell].shell;	if (shell[nshell].nnucl > 2)	  shell[nshell].nucl3 = i + shell[nshell].nucl3 - shell[nshell].shell;      }      shell[nshell].shell = i;      nshell++;    }  }  shfc->nshell = nshell;  shfc->shell  = shell;}

voiddo_redist_pos_coeffs(struct gmx_pme_t *pme, t_commrec *cr, int start, int homenr,                     gmx_bool bFirst, rvec x[], real *data){    int             d;    pme_atomcomm_t *atc;    atc = &pme->atc[0];    for (d = pme->ndecompdim - 1; d >= 0; d--)    {        int             n_d;        rvec           *x_d;        real           *param_d;        if (d == pme->ndecompdim - 1)        {            n_d     = homenr;            x_d     = x + start;            param_d = data;        }        else        {            n_d     = pme->atc[d + 1].n;            x_d     = atc->x;            param_d = atc->coefficient;        }        atc      = &pme->atc[d];        atc->npd = n_d;        if (atc->npd > atc->pd_nalloc)        {            atc->pd_nalloc = over_alloc_dd(atc->npd);            srenew(atc->pd, atc->pd_nalloc);        }        pme_calc_pidx_wrapper(n_d, pme->recipbox, x_d, atc);        where();        /* Redistribute x (only once) and qA/c6A or qB/c6B */        if (DOMAINDECOMP(cr))        {            dd_pmeredist_pos_coeffs(pme, n_d, bFirst, x_d, param_d, atc);        }    }}

real pull_potential(int ePull,t_pull *pull, t_mdatoms *md, t_pbc *pbc,		    t_commrec *cr, double t, real lambda,		    rvec *x, rvec *f, tensor vir, real *dvdlambda){  real V,dVdl;  pull_calc_coms(cr,pull,md,pbc,t,x,NULL);  do_pull_pot(ePull,pull,pbc,t,lambda,	      &V,pull->bVirial && MASTER(cr) ? vir : NULL,&dVdl);  /* Distribute forces over pulled groups */  apply_forces(pull, md, DOMAINDECOMP(cr) ? cr->dd->ga2la : NULL, f);  if (MASTER(cr)) {    *dvdlambda += dVdl;  }  return (MASTER(cr) ? V : 0.0);}

static void pull_set_pbcatom(t_commrec *cr, pull_group_work_t *pgrp,                             rvec *x,                             rvec x_pbc){    int a;    if (cr != NULL && DOMAINDECOMP(cr))    {        if (ga2la_get_home(cr->dd->ga2la, pgrp->params.pbcatom, &a))        {            copy_rvec(x[a], x_pbc);        }        else        {            clear_rvec(x_pbc);        }    }    else    {        copy_rvec(x[pgrp->params.pbcatom], x_pbc);    }}

//.........这里部分代码省略.........     * but we do so anyhow for consistency of the returned coordinates.     */    if (graph)    {        shift_self(graph, box, x);        if (TRICLINIC(box))        {            inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);        }        else        {            inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);        }    }    /* Check whether we need to do bondeds or correct for exclusions */    if (fr->bMolPBC &&        ((flags & GMX_FORCE_BONDED)         || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype)))    {        /* Since all atoms are in the rectangular or triclinic unit-cell,         * only single box vector shifts (2 in x) are required.         */        set_pbc_dd(&pbc, fr->ePBC, cr->dd, TRUE, box);    }    debug_gmx();    if (flags & GMX_FORCE_BONDED)    {        GMX_MPE_LOG(ev_calc_bonds_start);        wallcycle_sub_start(wcycle, ewcsBONDED);        calc_bonds(fplog, cr->ms,                   idef, x, hist, f, fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,                   DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,                   flags,                   fr->bSepDVDL && do_per_step(step, ir->nstlog), step);        /* Check if we have to determine energy differences         * at foreign lambda's.         */        if (fepvals->n_lambda > 0 && (flags & GMX_FORCE_DHDL) &&            idef->ilsort != ilsortNO_FE)        {            if (idef->ilsort != ilsortFE_SORTED)            {                gmx_incons("The bonded interactions are not sorted for free energy");            }            for (i = 0; i < enerd->n_lambda; i++)            {                reset_foreign_enerdata(enerd);                for (j = 0; j < efptNR; j++)                {                    lam_i[j] = (i == 0 ? lambda[j] : fepvals->all_lambda[j][i-1]);                }                calc_bonds_lambda(fplog, idef, x, fr, &pbc, graph, &(enerd->foreign_grpp), enerd->foreign_term, nrnb, lam_i, md,                                  fcd, DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);                sum_epot(&ir->opts, &(enerd->foreign_grpp), enerd->foreign_term);                enerd->enerpart_lambda[i] += enerd->foreign_term[F_EPOT];            }        }        debug_gmx();        GMX_MPE_LOG(ev_calc_bonds_finish);        wallcycle_sub_stop(wcycle, ewcsBONDED);    }    where();

//.........这里部分代码省略.........    /* Here sometimes we would not need to shift with NBFonly,     * but we do so anyhow for consistency of the returned coordinates.     */    if (graph)    {        shift_self(graph, box, x);        if (TRICLINIC(box))        {            inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);        }        else        {            inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);        }    }    /* Check whether we need to do listed interactions or correct for exclusions */    if (fr->bMolPBC &&        ((flags & GMX_FORCE_LISTED)         || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype)))    {        /* TODO There are no electrostatics methods that require this           transformation, when using the Verlet scheme, so update the           above conditional. */        /* Since all atoms are in the rectangular or triclinic unit-cell,         * only single box vector shifts (2 in x) are required.         */        set_pbc_dd(&pbc, fr->ePBC, cr->dd, TRUE, box);    }    debug_gmx();    do_force_listed(wcycle, box, ir->fepvals, cr->ms,                    idef, (const rvec *) x, hist, f, fr,                    &pbc, graph, enerd, nrnb, lambda, md, fcd,                    DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,                    flags);    where();    *cycles_pme = 0;    clear_mat(fr->vir_el_recip);    clear_mat(fr->vir_lj_recip);    /* Do long-range electrostatics and/or LJ-PME, including related short-range     * corrections.     */    if (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype))    {        int  status            = 0;        real Vlr_q             = 0, Vlr_lj = 0, Vcorr_q = 0, Vcorr_lj = 0;        real dvdl_long_range_q = 0, dvdl_long_range_lj = 0;        bSB = (ir->nwall == 2);        if (bSB)        {            copy_mat(box, boxs);            svmul(ir->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);            box_size[ZZ] *= ir->wall_ewald_zfac;        }        if (EEL_PME_EWALD(fr->eeltype) || EVDW_PME(fr->vdwtype))        {            real dvdl_long_range_correction_q   = 0;            real dvdl_long_range_correction_lj  = 0;            /* With the Verlet scheme exclusion forces are calculated             * in the non-bonded kernel.             */

void mdoutf_write_to_trajectory_files(FILE *fplog, t_commrec *cr,                                      gmx_mdoutf_t of,                                      int mdof_flags,                                      gmx_mtop_t *top_global,                                      gmx_int64_t step, double t,                                      t_state *state_local, t_state *state_global,                                      rvec *f_local, rvec *f_global){    rvec *local_v;    rvec *global_v;    /* MRS -- defining these variables is to manage the difference     * between half step and full step velocities, but there must be a better way . . . */    local_v  = state_local->v;    global_v = state_global->v;    if (DOMAINDECOMP(cr))    {        if (mdof_flags & MDOF_CPT)        {            dd_collect_state(cr->dd, state_local, state_global);        }        else        {            if (mdof_flags & (MDOF_X | MDOF_X_COMPRESSED))            {                dd_collect_vec(cr->dd, state_local, state_local->x,                               state_global->x);            }            if (mdof_flags & MDOF_V)            {                dd_collect_vec(cr->dd, state_local, local_v,                               global_v);            }        }        if (mdof_flags & MDOF_F)        {            dd_collect_vec(cr->dd, state_local, f_local, f_global);        }    }    else    {        if (mdof_flags & MDOF_CPT)        {            /* All pointers in state_local are equal to state_global,             * but we need to copy the non-pointer entries.             */            state_global->lambda = state_local->lambda;            state_global->veta   = state_local->veta;            state_global->vol0   = state_local->vol0;            copy_mat(state_local->box, state_global->box);            copy_mat(state_local->boxv, state_global->boxv);            copy_mat(state_local->svir_prev, state_global->svir_prev);            copy_mat(state_local->fvir_prev, state_global->fvir_prev);            copy_mat(state_local->pres_prev, state_global->pres_prev);        }    }    if (MASTER(cr))    {        if (mdof_flags & MDOF_CPT)        {            fflush_tng(of->tng);            fflush_tng(of->tng_low_prec);            write_checkpoint(of->fn_cpt, of->bKeepAndNumCPT,                             fplog, cr, of->eIntegrator, of->simulation_part,                             of->bExpanded, of->elamstats, step, t, state_global);        }        if (mdof_flags & (MDOF_X | MDOF_V | MDOF_F))        {            if (of->fp_trn)            {                gmx_trr_write_frame(of->fp_trn, step, t, state_local->lambda[efptFEP],                                    state_local->box, top_global->natoms,                                    (mdof_flags & MDOF_X) ? state_global->x : NULL,                                    (mdof_flags & MDOF_V) ? global_v : NULL,                                    (mdof_flags & MDOF_F) ? f_global : NULL);                if (gmx_fio_flush(of->fp_trn) != 0)                {                    gmx_file("Cannot write trajectory; maybe you are out of disk space?");                }            }            gmx_fwrite_tng(of->tng, FALSE, step, t, state_local->lambda[efptFEP],                           state_local->box,                           top_global->natoms,                           (mdof_flags & MDOF_X) ? state_global->x : NULL,                           (mdof_flags & MDOF_V) ? global_v : NULL,                           (mdof_flags & MDOF_F) ? f_global : NULL);        }        if (mdof_flags & MDOF_X_COMPRESSED)        {            rvec *xxtc = NULL;            if (of->natoms_x_compressed == of->natoms_global)            {                /* We are writing the positions of all of the atoms to                   the compressed output *///.........这里部分代码省略.........

void pme_loadbal_do(pme_load_balancing_t *pme_lb,                    t_commrec            *cr,                    FILE                 *fp_err,                    FILE                 *fp_log,                    t_inputrec           *ir,                    t_forcerec           *fr,                    t_state              *state,                    gmx_wallcycle_t       wcycle,                    gmx_int64_t           step,                    gmx_int64_t           step_rel,                    gmx_bool             *bPrinting){    int    n_prev;    double cycles_prev;    assert(pme_lb != NULL);    if (!pme_lb->bActive)    {        return;    }    n_prev      = pme_lb->cycles_n;    cycles_prev = pme_lb->cycles_c;    wallcycle_get(wcycle, ewcSTEP, &pme_lb->cycles_n, &pme_lb->cycles_c);    if (pme_lb->cycles_n == 0)    {        /* Before the first step we haven't done any steps yet */        return;    }    /* Sanity check, we expect nstlist cycle counts */    if (pme_lb->cycles_n - n_prev != ir->nstlist)    {        /* We could return here, but it's safer to issue and error and quit */        gmx_incons("pme_loadbal_do called at an interval != nstlist");    }    /* PME grid + cut-off optimization with GPUs or PME ranks */    if (!pme_lb->bBalance && pme_lb->bSepPMERanks)    {        if (pme_lb->bTriggerOnDLB)        {            pme_lb->bBalance = dd_dlb_is_on(cr->dd);        }        /* We should ignore the first timing to avoid timing allocation         * overhead. And since the PME load balancing is called just         * before DD repartitioning, the ratio returned by dd_pme_f_ratio         * is not over the last nstlist steps, but the nstlist steps before         * that. So the first useful ratio is available at step_rel=3*nstlist.         */        else if (step_rel >= 3*ir->nstlist)        {            if (DDMASTER(cr->dd))            {                /* If PME rank load is too high, start tuning */                pme_lb->bBalance =                    (dd_pme_f_ratio(cr->dd) >= loadBalanceTriggerFactor);            }            dd_bcast(cr->dd, sizeof(gmx_bool), &pme_lb->bBalance);        }        pme_lb->bActive = (pme_lb->bBalance ||                           step_rel <= pme_lb->step_rel_stop);    }    /* The location in the code of this balancing termination is strange.     * You would expect to have it after the call to pme_load_balance()     * below, since there pme_lb->stage is updated.     * But when terminating directly after deciding on and selecting the     * optimal setup, DLB will turn on right away if it was locked before.     * This might be due to PME reinitialization. So we check stage here     * to allow for another nstlist steps with DLB locked to stabilize     * the performance.     */    if (pme_lb->bBalance && pme_lb->stage == pme_lb->nstage)    {        pme_lb->bBalance = FALSE;        if (DOMAINDECOMP(cr) && dd_dlb_is_locked(cr->dd))        {            /* Unlock the DLB=auto, DLB is allowed to activate */            dd_dlb_unlock(cr->dd);            md_print_warn(cr, fp_log, "NOTE: DLB can now turn on, when beneficial/n");            /* We don't deactivate the tuning yet, since we will balance again             * after DLB gets turned on, if it does within PMETune_period.             */            continue_pme_loadbal(pme_lb, TRUE);            pme_lb->bTriggerOnDLB = TRUE;            pme_lb->step_rel_stop = step_rel + PMETunePeriod*ir->nstlist;        }        else        {            /* We're completely done with PME tuning */            pme_lb->bActive = FALSE;        }        if (DOMAINDECOMP(cr))        {            /* Set the cut-off limit to the final selected cut-off,//.........这里部分代码省略.........

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

static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,                             t_pbc *pbc, double t, rvec *x, rvec *xp){    int           c, i, ii, m, start, end;    rvec          g_x, dx, dir;    double        r0_2, sum_a, sum_ap, dr2, mass, weight, wmass, wwmass, inp;    t_pull_coord *pcrd;    t_pull_group *pref, *pgrp, *pdyna;    gmx_ga2la_t   ga2la = NULL;    if (pull->dbuf_cyl == NULL)    {        snew(pull->dbuf_cyl, pull->ncoord*4);    }    if (cr && DOMAINDECOMP(cr))    {        ga2la = cr->dd->ga2la;    }    start = 0;    end   = md->homenr;    r0_2 = dsqr(pull->cyl_r0);    /* loop over all groups to make a reference group for each*/    for (c = 0; c < pull->ncoord; c++)    {        pcrd  = &pull->coord[c];        /* pref will be the same group for all pull coordinates */        pref  = &pull->group[pcrd->group[0]];        pgrp  = &pull->group[pcrd->group[1]];        pdyna = &pull->dyna[c];        copy_rvec(pcrd->vec, dir);        sum_a          = 0;        sum_ap         = 0;        wmass          = 0;        wwmass         = 0;        pdyna->nat_loc = 0;        for (m = 0; m < DIM; m++)        {            g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);        }        /* loop over all atoms in the main ref group */        for (i = 0; i < pref->nat; i++)        {            ii = pref->ind[i];            if (ga2la)            {                if (!ga2la_get_home(ga2la, pref->ind[i], &ii))                {                    ii = -1;                }            }            if (ii >= start && ii < end)            {                pbc_dx_aiuc(pbc, x[ii], g_x, dx);                inp = iprod(dir, dx);                dr2 = 0;                for (m = 0; m < DIM; m++)                {                    dr2 += dsqr(dx[m] - inp*dir[m]);                }                if (dr2 < r0_2)                {                    /* add to index, to sum of COM, to weight array */                    if (pdyna->nat_loc >= pdyna->nalloc_loc)                    {                        pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);                        srenew(pdyna->ind_loc, pdyna->nalloc_loc);                        srenew(pdyna->weight_loc, pdyna->nalloc_loc);                    }                    pdyna->ind_loc[pdyna->nat_loc] = ii;                    mass   = md->massT[ii];                    weight = get_weight(sqrt(dr2), pull->cyl_r1, pull->cyl_r0);                    pdyna->weight_loc[pdyna->nat_loc] = weight;                    sum_a += mass*weight*inp;                    if (xp)                    {                        pbc_dx_aiuc(pbc, xp[ii], g_x, dx);                        inp     = iprod(dir, dx);                        sum_ap += mass*weight*inp;                    }                    wmass  += mass*weight;                    wwmass += mass*sqr(weight);                    pdyna->nat_loc++;                }            }        }        pull->dbuf_cyl[c*4+0] = wmass;        pull->dbuf_cyl[c*4+1] = wwmass;        pull->dbuf_cyl[c*4+2] = sum_a;        pull->dbuf_cyl[c*4+3] = sum_ap;    }    if (cr && PAR(cr))//.........这里部分代码省略.........

static void do_lincs(rvec *x,rvec *xp,matrix box,t_pbc *pbc,                     struct gmx_lincsdata *lincsd,real *invmass,					 t_commrec *cr,                     real wangle,int *warn,                     real invdt,rvec *v,                     gmx_bool bCalcVir,tensor rmdr){    int     b,i,j,k,n,iter;    real    tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,len2,dlen2,wfac,lam;      rvec    dx;    int     ncons,*bla,*blnr,*blbnb;    rvec    *r;    real    *blc,*blmf,*bllen,*blcc,*rhs1,*rhs2,*sol,*lambda;    int     *nlocat;        ncons  = lincsd->nc;    bla    = lincsd->bla;    r      = lincsd->tmpv;    blnr   = lincsd->blnr;    blbnb  = lincsd->blbnb;    blc    = lincsd->blc;    blmf   = lincsd->blmf;    bllen  = lincsd->bllen;    blcc   = lincsd->tmpncc;    rhs1   = lincsd->tmp1;    rhs2   = lincsd->tmp2;    sol    = lincsd->tmp3;    lambda = lincsd->lambda;        if (DOMAINDECOMP(cr) && cr->dd->constraints)    {        nlocat = dd_constraints_nlocalatoms(cr->dd);    }    else if (PARTDECOMP(cr))    {        nlocat = pd_constraints_nlocalatoms(cr->pd);    }    else    {        nlocat = NULL;    }        *warn = 0;    if (pbc)    {        /* Compute normalized i-j vectors */        for(b=0; b<ncons; b++)        {            pbc_dx_aiuc(pbc,x[bla[2*b]],x[bla[2*b+1]],dx);            unitv(dx,r[b]);        }          for(b=0; b<ncons; b++)        {            for(n=blnr[b]; n<blnr[b+1]; n++)            {                blcc[n] = blmf[n]*iprod(r[b],r[blbnb[n]]);            }            pbc_dx_aiuc(pbc,xp[bla[2*b]],xp[bla[2*b+1]],dx);            mvb = blc[b]*(iprod(r[b],dx) - bllen[b]);            rhs1[b] = mvb;            sol[b]  = mvb;        }    }    else    {        /* Compute normalized i-j vectors */        for(b=0; b<ncons; b++)        {            i = bla[2*b];            j = bla[2*b+1];            tmp0 = x[i][0] - x[j][0];            tmp1 = x[i][1] - x[j][1];            tmp2 = x[i][2] - x[j][2];            rlen = gmx_invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);            r[b][0] = rlen*tmp0;            r[b][1] = rlen*tmp1;            r[b][2] = rlen*tmp2;        } /* 16 ncons flops */                for(b=0; b<ncons; b++)        {            tmp0 = r[b][0];            tmp1 = r[b][1];            tmp2 = r[b][2];            len = bllen[b];            i = bla[2*b];            j = bla[2*b+1];            for(n=blnr[b]; n<blnr[b+1]; n++)            {                k = blbnb[n];                blcc[n] = blmf[n]*(tmp0*r[k][0] + tmp1*r[k][1] + tmp2*r[k][2]);             } /* 6 nr flops */            mvb = blc[b]*(tmp0*(xp[i][0] - xp[j][0]) +                          tmp1*(xp[i][1] - xp[j][1]) +                              tmp2*(xp[i][2] - xp[j][2]) - len);            rhs1[b] = mvb;            sol[b]  = mvb;            /* 10 flops */        }//.........这里部分代码省略.........

gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,                          t_state *state,real *ener,                          t_state *state_local,                          int step,real time){    gmx_multisim_t *ms;    int  exchange=-1,shift;    gmx_bool bExchanged=FALSE;        ms = cr->ms;      if (MASTER(cr))    {        exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box),                                        step,time);        bExchanged = (exchange >= 0);    }        if (PAR(cr))    {#ifdef GMX_MPI        MPI_Bcast(&bExchanged,sizeof(gmx_bool),MPI_BYTE,MASTERRANK(cr),                  cr->mpi_comm_mygroup);#endif    }        if (bExchanged)    {        /* Exchange the states */        if (PAR(cr))        {            /* Collect the global state on the master node */            if (DOMAINDECOMP(cr))            {                dd_collect_state(cr->dd,state_local,state);            }            else            {                pd_collect_state(cr,state);            }        }                if (MASTER(cr))        {            /* Exchange the global states between the master nodes */            if (debug)            {                fprintf(debug,"Exchanging %d with %d/n",ms->sim,exchange);            }            exchange_state(ms,exchange,state);                        if (re->type == ereTEMP)            {                scale_velocities(state,sqrt(re->q[ms->sim]/re->q[exchange]));            }        }        /* With domain decomposition the global state is distributed later */        if (!DOMAINDECOMP(cr))        {            /* Copy the global state to the local state data structure */            copy_state_nonatomdata(state,state_local);                        if (PAR(cr))            {                bcast_state(cr,state,FALSE);            }        }    }            return bExchanged;}

/* calculates center of mass of selection index from all coordinates x */void pull_calc_coms(t_commrec *cr,                    t_pull *pull, t_mdatoms *md, t_pbc *pbc, double t,                    rvec x[], rvec *xp){    int           g, i, ii, m;    real          mass, w, wm, twopi_box = 0;    double        wmass, wwmass, invwmass;    dvec          com, comp;    double        cm, sm, cmp, smp, ccm, csm, ssm, csw, snw;    rvec         *xx[2], x_pbc = {0, 0, 0}, dx;    t_pull_group *pgrp;    if (pull->rbuf == NULL)    {        snew(pull->rbuf, pull->ngroup);    }    if (pull->dbuf == NULL)    {        snew(pull->dbuf, 3*pull->ngroup);    }    if (pull->bRefAt && pull->bSetPBCatoms)    {        pull_set_pbcatoms(cr, pull, x, pull->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)    {        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++)    {        pgrp = &pull->group[g];        clear_dvec(com);        clear_dvec(comp);        wmass  = 0;        wwmass = 0;        cm     = 0;        sm     = 0;        cmp    = 0;        smp    = 0;        ccm    = 0;        csm    = 0;        ssm    = 0;        if (!(g == 0 && PULL_CYL(pull)))        {            if (pgrp->epgrppbc == epgrppbcREFAT)            {                /* Set the pbc atom */                copy_rvec(pull->rbuf[g], x_pbc);            }            w = 1;            for (i = 0; i < pgrp->nat_loc; i++)            {                ii   = pgrp->ind_loc[i];                mass = md->massT[ii];                if (pgrp->epgrppbc != epgrppbcCOS)                {                    if (pgrp->weight_loc)                    {                        w = pgrp->weight_loc[i];                    }                    wm      = w*mass;                    wmass  += wm;                    wwmass += wm*w;                    if (pgrp->epgrppbc == epgrppbcNONE)                    {                        /* Plain COM: sum the coordinates */                        for (m = 0; m < DIM; m++)                        {                            com[m]    += wm*x[ii][m];                        }                        if (xp)                        {                            for (m = 0; m < DIM; m++)                            {                                comp[m] += wm*xp[ii][m];                            }                        }//.........这里部分代码省略.........

//.........这里部分代码省略.........    /* Here sometimes we would not need to shift with NBFonly,     * but we do so anyhow for consistency of the returned coordinates.     */    if (graph)    {        shift_self(graph,box,x);        if (TRICLINIC(box))        {            inc_nrnb(nrnb,eNR_SHIFTX,2*graph->nnodes);        }        else        {            inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);        }    }    /* Check whether we need to do bondeds or correct for exclusions */    if (fr->bMolPBC &&            ((flags & GMX_FORCE_BONDED)             || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype)))    {        /* Since all atoms are in the rectangular or triclinic unit-cell,         * only single box vector shifts (2 in x) are required.         */        set_pbc_dd(&pbc,fr->ePBC,cr->dd,TRUE,box);    }    debug_gmx();    if (flags & GMX_FORCE_BONDED)    {        GMX_MPE_LOG(ev_calc_bonds_start);        calc_bonds(fplog,cr->ms,                   idef,x,hist,f,fr,&pbc,graph,enerd,nrnb,lambda,md,fcd,                   DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL, atype, born,                   fr->bSepDVDL && do_per_step(step,ir->nstlog),step);        /* Check if we have to determine energy differences         * at foreign lambda's.         */        if (ir->n_flambda > 0 && (flags & GMX_FORCE_DHDL) &&                idef->ilsort != ilsortNO_FE)        {            if (idef->ilsort != ilsortFE_SORTED)            {                gmx_incons("The bonded interactions are not sorted for free energy");            }            init_enerdata(mtop->groups.grps[egcENER].nr,ir->n_flambda,&ed_lam);            for(i=0; i<enerd->n_lambda; i++)            {                lam_i = (i==0 ? lambda : ir->flambda[i-1]);                dvdl_dum = 0;                reset_enerdata(&ir->opts,fr,TRUE,&ed_lam,FALSE);                calc_bonds_lambda(fplog,                                  idef,x,fr,&pbc,graph,&ed_lam,nrnb,lam_i,md,                                  fcd,                                  DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);                sum_epot(&ir->opts,&ed_lam);                enerd->enerpart_lambda[i] += ed_lam.term[F_EPOT];            }            destroy_enerdata(&ed_lam);        }        debug_gmx();        GMX_MPE_LOG(ev_calc_bonds_finish);    }

//.........这里部分代码省略.........        for (j = 0; (j < egNR); j++)        {            inn[j] = add_binr(rb, enerd->grpp.nener, enerd->grpp.ener[j]);        }        where();        if (inputrec->efep != efepNO)        {            idvdll  = add_bind(rb, efptNR, enerd->dvdl_lin);            idvdlnl = add_bind(rb, efptNR, enerd->dvdl_nonlin);            if (enerd->n_lambda > 0)            {                iepl = add_bind(rb, enerd->n_lambda, enerd->enerpart_lambda);            }        }    }    if (vcm)    {        icm   = add_binr(rb, DIM*vcm->nr, vcm->group_p[0]);        where();        imass = add_binr(rb, vcm->nr, vcm->group_mass);        where();        if (vcm->mode == ecmANGULAR)        {            icj   = add_binr(rb, DIM*vcm->nr, vcm->group_j[0]);            where();            icx   = add_binr(rb, DIM*vcm->nr, vcm->group_x[0]);            where();            ici   = add_binr(rb, DIM*DIM*vcm->nr, vcm->group_i[0][0]);            where();        }    }    if (DOMAINDECOMP(cr))    {        nb  = cr->dd->nbonded_local;        inb = add_bind(rb, 1, &nb);    }    where();    if (nsig > 0)    {        isig = add_binr(rb, nsig, sig);    }    /* Global sum it all */    if (debug)    {        fprintf(debug, "Summing %d energies/n", rb->maxreal);    }    sum_bin(rb, cr);    where();    /* Extract all the data locally */    if (bConstrVir)    {        extract_binr(rb, isv, DIM*DIM, svir[0]);    }    /* We need the force virial and the kinetic energy for the first time through with velocity verlet */    if (bTemp || !bVV)    {        if (ekind)        {            for (j = 0; (j < inputrec->opts.ngtc); j++)            {

gmx_bool pme_load_balance(pme_load_balancing_t       pme_lb,                          t_commrec                 *cr,                          FILE                      *fp_err,                          FILE                      *fp_log,                          t_inputrec                *ir,                          t_state                   *state,                          double                     cycles,                          interaction_const_t       *ic,                          struct nonbonded_verlet_t *nbv,                          struct gmx_pme_t **        pmedata,                          gmx_int64_t                step){    gmx_bool     OK;    pme_setup_t *set;    double       cycles_fast;    char         buf[STRLEN], sbuf[22];    real         rtab;    gmx_bool     bUsesSimpleTables = TRUE;    if (pme_lb->stage == pme_lb->nstage)    {        return FALSE;    }    if (PAR(cr))    {        gmx_sumd(1, &cycles, cr);        cycles /= cr->nnodes;    }    set = &pme_lb->setup[pme_lb->cur];    set->count++;    rtab = ir->rlistlong + ir->tabext;    if (set->count % 2 == 1)    {        /* Skip the first cycle, because the first step after a switch         * is much slower due to allocation and/or caching effects.         */        return TRUE;    }    sprintf(buf, "step %4s: ", gmx_step_str(step, sbuf));    print_grid(fp_err, fp_log, buf, "timed with", set, cycles);    if (set->count <= 2)    {        set->cycles = cycles;    }    else    {        if (cycles*PME_LB_ACCEL_TOL < set->cycles &&            pme_lb->stage == pme_lb->nstage - 1)        {            /* The performance went up a lot (due to e.g. DD load balancing).             * Add a stage, keep the minima, but rescan all setups.             */            pme_lb->nstage++;            if (debug)            {                fprintf(debug, "The performance for grid %d %d %d went from %.3f to %.1f M-cycles, this is more than %f/n"                        "Increased the number stages to %d"                        " and ignoring the previous performance/n",                        set->grid[XX], set->grid[YY], set->grid[ZZ],                        cycles*1e-6, set->cycles*1e-6, PME_LB_ACCEL_TOL,                        pme_lb->nstage);            }        }        set->cycles = min(set->cycles, cycles);    }    if (set->cycles < pme_lb->setup[pme_lb->fastest].cycles)    {        pme_lb->fastest = pme_lb->cur;        if (DOMAINDECOMP(cr))        {            /* We found a new fastest setting, ensure that with subsequent             * shorter cut-off's the dynamic load balancing does not make             * the use of the current cut-off impossible. This solution is             * a trade-off, as the PME load balancing and DD domain size             * load balancing can interact in complex ways.             * With the Verlet kernels, DD load imbalance will usually be             * mainly due to bonded interaction imbalance, which will often             * quickly push the domain boundaries beyond the limit for the             * optimal, PME load balanced, cut-off. But it could be that             * better overal performance can be obtained with a slightly             * shorter cut-off and better DD load balancing.             */            change_dd_dlb_cutoff_limit(cr);        }    }    cycles_fast = pme_lb->setup[pme_lb->fastest].cycles;    /* Check in stage 0 if we should stop scanning grids.     * Stop when the time is more than SLOW_FAC longer than the fastest.     */    if (pme_lb->stage == 0 && pme_lb->cur > 0 &&//.........这里部分代码省略.........

void update_QMMMrec(t_commrec *cr,		    t_forcerec *fr,		    rvec x[],		    t_mdatoms *md,		    matrix box,		    gmx_localtop_t *top){  /* updates the coordinates of both QM atoms and MM atoms and stores   * them in the QMMMrec.   *   * NOTE: is NOT yet working if there are no PBC. Also in ns.c, simple   * ns needs to be fixed!   */  int    mm_max=0,mm_nr=0,mm_nr_new,i,j,is,k,shift;  t_j_particle    *mm_j_particles=NULL,*qm_i_particles=NULL;  t_QMMMrec    *qr;  t_nblist    QMMMlist;  rvec    dx,crd;  int    *MMatoms;  t_QMrec    *qm;  t_MMrec    *mm;  t_pbc    pbc;  int    *parallelMMarray=NULL;  real    c12au,c6au;  c6au  = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM,6));  c12au = (HARTREE2KJ*AVOGADRO*pow(BOHR2NM,12));  /* every cpu has this array. On every processor we fill this array   * with 1's and 0's. 1's indicate the atoms is a QM atom on the   * current cpu in a later stage these arrays are all summed. indexes   * > 0 indicate the atom is a QM atom. Every node therefore knows   * whcih atoms are part of the QM subsystem.   */  /* copy some pointers */  qr          = fr->qr;  mm          = qr->mm;  QMMMlist    = fr->QMMMlist;  /*  init_pbc(box);  needs to be called first, see pbc.h */  set_pbc_dd(&pbc,fr->ePBC,DOMAINDECOMP(cr) ? cr->dd : NULL,FALSE,box);  /* only in standard (normal) QMMM we need the neighbouring MM   * particles to provide a electric field of point charges for the QM   * atoms.   */  if(qr->QMMMscheme==eQMMMschemenormal){ /* also implies 1 QM-layer */    /* we NOW create/update a number of QMMMrec entries:     *     * 1) the shiftQM, containing the shifts of the QM atoms     *     * 2) the indexMM array, containing the index of the MM atoms     *     * 3) the shiftMM, containing the shifts of the MM atoms     *     * 4) the shifted coordinates of the MM atoms     *     * the shifts are used for computing virial of the QM/MM particles.     */    qm = qr->qm[0]; /* in case of normal QMMM, there is only one group */    snew(qm_i_particles,QMMMlist.nri);    if(QMMMlist.nri){      qm_i_particles[0].shift = XYZ2IS(0,0,0);      for(i=0;i<QMMMlist.nri;i++){	qm_i_particles[i].j     = QMMMlist.iinr[i];	if(i){	  qm_i_particles[i].shift = pbc_dx_aiuc(&pbc,x[QMMMlist.iinr[0]],						x[QMMMlist.iinr[i]],dx);	}	/* However, since nri >= nrQMatoms, we do a quicksort, and throw	 * out double, triple, etc. entries later, as we do for the MM	 * list too.	 */	/* compute the shift for the MM j-particles with respect to	 * the QM i-particle and store them.	 */	crd[0] = IS2X(QMMMlist.shift[i]) + IS2X(qm_i_particles[i].shift);	crd[1] = IS2Y(QMMMlist.shift[i]) + IS2Y(qm_i_particles[i].shift);	crd[2] = IS2Z(QMMMlist.shift[i]) + IS2Z(qm_i_particles[i].shift);	is = XYZ2IS(crd[0],crd[1],crd[2]);	for(j=QMMMlist.jindex[i];	    j<QMMMlist.jindex[i+1];	    j++){	  if(mm_nr >= mm_max){//.........这里部分代码省略.........

void check_resource_division_efficiency(const gmx_hw_info_t *hwinfo,                                        const gmx_hw_opt_t  *hw_opt,                                        gmx_bool             bNtOmpOptionSet,                                        t_commrec           *cr,                                        FILE                *fplog){#if defined GMX_OPENMP && defined GMX_MPI    int         nth_omp_min, nth_omp_max, ngpu;    char        buf[1000];#ifdef GMX_THREAD_MPI    const char *mpi_option = " (option -ntmpi)";#else    const char *mpi_option = "";#endif    /* This function should be called after thread-MPI (when configured) and     * OpenMP have been initialized. Check that here.     */#ifdef GMX_THREAD_MPI    GMX_RELEASE_ASSERT(nthreads_omp_faster_default >= nthreads_omp_mpi_ok_max, "Inconsistent OpenMP thread count default values");    GMX_RELEASE_ASSERT(hw_opt->nthreads_tmpi >= 1, "Must have at least one thread-MPI rank");#endif    GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault) >= 1, "Must have at least one OpenMP thread");    nth_omp_min = gmx_omp_nthreads_get(emntDefault);    nth_omp_max = gmx_omp_nthreads_get(emntDefault);    ngpu        = hw_opt->gpu_opt.n_dev_use;    /* Thread-MPI seems to have a bug with reduce on 1 node, so use a cond. */    if (cr->nnodes + cr->npmenodes > 1)    {        int count[3], count_max[3];        count[0] = -nth_omp_min;        count[1] =  nth_omp_max;        count[2] =  ngpu;        MPI_Allreduce(count, count_max, 3, MPI_INT, MPI_MAX, cr->mpi_comm_mysim);        /* In case of an inhomogeneous run setup we use the maximum counts */        nth_omp_min = -count_max[0];        nth_omp_max =  count_max[1];        ngpu        =  count_max[2];    }    int nthreads_omp_mpi_ok_min;    if (ngpu == 0)    {        nthreads_omp_mpi_ok_min = nthreads_omp_mpi_ok_min_cpu;    }    else    {        /* With GPUs we set the minimum number of OpenMP threads to 2 to catch         * cases where the user specifies #ranks == #cores.         */        nthreads_omp_mpi_ok_min = nthreads_omp_mpi_ok_min_gpu;    }    if (DOMAINDECOMP(cr) && cr->nnodes > 1)    {        if (nth_omp_max < nthreads_omp_mpi_ok_min ||            (!(ngpu > 0 && !gmx_gpu_sharing_supported()) &&             nth_omp_max > nthreads_omp_mpi_ok_max))        {            /* Note that we print target_max here, not ok_max */            sprintf(buf, "Your choice of number of MPI ranks and amount of resources results in using %d OpenMP threads per rank, which is most likely inefficient. The optimum is usually between %d and %d threads per rank.",                    nth_omp_max,                    nthreads_omp_mpi_ok_min,                    nthreads_omp_mpi_target_max);            if (bNtOmpOptionSet)            {                md_print_warn(cr, fplog, "NOTE: %s/n", buf);            }            else            {                /* This fatal error, and the one below, is nasty, but it's                 * probably the only way to ensure that all users don't waste                 * a lot of resources, since many users don't read logs/stderr.                 */                gmx_fatal(FARGS, "%s If you want to run with this setup, specify the -ntomp option. But we suggest to change the number of MPI ranks%s.", buf, mpi_option);            }        }    }    else    {        /* No domain decomposition (or only one domain) */        if (!(ngpu > 0 && !gmx_gpu_sharing_supported()) &&            nth_omp_max > nthreads_omp_faster(hwinfo->cpuid_info, ngpu > 0))        {            /* To arrive here, the user/system set #ranks and/or #OMPthreads */            gmx_bool bEnvSet;            char     buf2[256];            bEnvSet = (getenv("OMP_NUM_THREADS") != NULL);            if (bNtOmpOptionSet || bEnvSet)            {                sprintf(buf2, "You requested %d OpenMP threads", nth_omp_max);//.........这里部分代码省略.........

//.........这里部分代码省略.........    /* Here sometimes we would not need to shift with NBFonly,     * but we do so anyhow for consistency of the returned coordinates.     */    if (graph)    {        shift_self(graph, box, x);        if (TRICLINIC(box))        {            inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);        }        else        {            inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);        }    }    /* Check whether we need to do listed interactions or correct for exclusions */    if (fr->bMolPBC &&        ((flags & GMX_FORCE_LISTED)         || EEL_RF(fr->eeltype) || EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype)))    {        /* TODO There are no electrostatics methods that require this           transformation, when using the Verlet scheme, so update the           above conditional. */        /* Since all atoms are in the rectangular or triclinic unit-cell,         * only single box vector shifts (2 in x) are required.         */        set_pbc_dd(&pbc, fr->ePBC, cr->dd, TRUE, box);    }    debug_gmx();    do_force_listed(wcycle, box, ir->fepvals, cr->ms,                    idef, (const rvec *) x, hist, f, fr,                    &pbc, graph, enerd, nrnb, lambda, md, fcd,                    DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL,                    flags);    where();    *cycles_pme = 0;    clear_mat(fr->vir_el_recip);    clear_mat(fr->vir_lj_recip);    /* Do long-range electrostatics and/or LJ-PME, including related short-range     * corrections.     */    if (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype))    {        int  status            = 0;        real Vlr_q             = 0, Vlr_lj = 0, Vcorr_q = 0, Vcorr_lj = 0;        real dvdl_long_range_q = 0, dvdl_long_range_lj = 0;        bSB = (ir->nwall == 2);        if (bSB)        {            copy_mat(box, boxs);            svmul(ir->wall_ewald_zfac, boxs[ZZ], boxs[ZZ]);            box_size[ZZ] *= ir->wall_ewald_zfac;        }        if (EEL_PME_EWALD(fr->eeltype) || EVDW_PME(fr->vdwtype))        {            real dvdl_long_range_correction_q   = 0;            real dvdl_long_range_correction_lj  = 0;            /* With the Verlet scheme exclusion forces are calculated             * in the non-bonded kernel.             */

gmx_bool pme_load_balance(pme_load_balancing_t pme_lb,                          t_commrec           *cr,                          FILE                *fp_err,                          FILE                *fp_log,                          t_inputrec          *ir,                          t_state             *state,                          double               cycles,                          interaction_const_t *ic,                          nonbonded_verlet_t  *nbv,                          gmx_pme_t           *pmedata,                          gmx_large_int_t      step){    gmx_bool     OK;    pme_setup_t *set;    double       cycles_fast;    char         buf[STRLEN], sbuf[22];    real         rtab;    gmx_bool     bUsesSimpleTables = TRUE;    if (pme_lb->stage == pme_lb->nstage)    {        return FALSE;    }    if (PAR(cr))    {        gmx_sumd(1, &cycles, cr);        cycles /= cr->nnodes;    }    set = &pme_lb->setup[pme_lb->cur];    set->count++;    rtab = ir->rlistlong + ir->tabext;    if (set->count % 2 == 1)    {        /* Skip the first cycle, because the first step after a switch         * is much slower due to allocation and/or caching effects.         */        return TRUE;    }    sprintf(buf, "step %4s: ", gmx_step_str(step, sbuf));    print_grid(fp_err, fp_log, buf, "timed with", set, cycles);    if (set->count <= 2)    {        set->cycles = cycles;    }    else    {        if (cycles*PME_LB_ACCEL_TOL < set->cycles &&            pme_lb->stage == pme_lb->nstage - 1)        {            /* The performance went up a lot (due to e.g. DD load balancing).             * Add a stage, keep the minima, but rescan all setups.             */            pme_lb->nstage++;            if (debug)            {                fprintf(debug, "The performance for grid %d %d %d went from %.3f to %.1f M-cycles, this is more than %f/n"                        "Increased the number stages to %d"                        " and ignoring the previous performance/n",                        set->grid[XX], set->grid[YY], set->grid[ZZ],                        cycles*1e-6, set->cycles*1e-6, PME_LB_ACCEL_TOL,                        pme_lb->nstage);            }        }        set->cycles = min(set->cycles, cycles);    }    if (set->cycles < pme_lb->setup[pme_lb->fastest].cycles)    {        pme_lb->fastest = pme_lb->cur;        if (DOMAINDECOMP(cr))        {            /* We found a new fastest setting, ensure that with subsequent             * shorter cut-off's the dynamic load balancing does not make             * the use of the current cut-off impossible. This solution is             * a trade-off, as the PME load balancing and DD domain size             * load balancing can interact in complex ways.             * With the Verlet kernels, DD load imbalance will usually be             * mainly due to bonded interaction imbalance, which will often             * quickly push the domain boundaries beyond the limit for the             * optimal, PME load balanced, cut-off. But it could be that             * better overal performance can be obtained with a slightly             * shorter cut-off and better DD load balancing.             */            change_dd_dlb_cutoff_limit(cr);        }    }    cycles_fast = pme_lb->setup[pme_lb->fastest].cycles;    /* Check in stage 0 if we should stop scanning grids.     * Stop when the time is more than SLOW_FAC longer than the fastest.     */    if (pme_lb->stage == 0 && pme_lb->cur > 0 &&//.........这里部分代码省略.........

void set_lincs(t_idef *idef,t_mdatoms *md,               gmx_bool bDynamics,t_commrec *cr,               struct gmx_lincsdata *li){    int      start,natoms,nflexcon;    t_blocka at2con;    t_iatom  *iatom;    int      i,k,ncc_alloc,ni,con,nconnect,concon;    int      type,a1,a2;    real     lenA=0,lenB;    gmx_bool     bLocal;    li->nc = 0;    li->ncc = 0;		    /* This is the local topology, so there are only F_CONSTR constraints */    if (idef->il[F_CONSTR].nr == 0)    {        /* There are no constraints,         * we do not need to fill any data structures.         */        return;    }        if (debug)    {        fprintf(debug,"Building the LINCS connectivity/n");    }        if (DOMAINDECOMP(cr))    {        if (cr->dd->constraints)        {            dd_get_constraint_range(cr->dd,&start,&natoms);        }        else        {            natoms = cr->dd->nat_home;        }        start = 0;    }    else if(PARTDECOMP(cr))	{		pd_get_constraint_range(cr->pd,&start,&natoms);	}	else    {        start  = md->start;        natoms = md->homenr;    }    at2con = make_at2con(start,natoms,idef->il,idef->iparams,bDynamics,                         &nflexcon);	    if (idef->il[F_CONSTR].nr/3 > li->nc_alloc || li->nc_alloc == 0)    {        li->nc_alloc = over_alloc_dd(idef->il[F_CONSTR].nr/3);        srenew(li->bllen0,li->nc_alloc);        srenew(li->ddist,li->nc_alloc);        srenew(li->bla,2*li->nc_alloc);        srenew(li->blc,li->nc_alloc);        srenew(li->blc1,li->nc_alloc);        srenew(li->blnr,li->nc_alloc+1);        srenew(li->bllen,li->nc_alloc);        srenew(li->tmpv,li->nc_alloc);        srenew(li->tmp1,li->nc_alloc);        srenew(li->tmp2,li->nc_alloc);        srenew(li->tmp3,li->nc_alloc);        srenew(li->lambda,li->nc_alloc);        if (li->ncg_triangle > 0)        {            /* This is allocating too much, but it is difficult to improve */            srenew(li->triangle,li->nc_alloc);            srenew(li->tri_bits,li->nc_alloc);        }    }        iatom = idef->il[F_CONSTR].iatoms;        ncc_alloc = li->ncc_alloc;    li->blnr[0] = 0;        ni = idef->il[F_CONSTR].nr/3;    con = 0;    nconnect = 0;    li->blnr[con] = nconnect;    for(i=0; i<ni; i++)    {        bLocal = TRUE;        type = iatom[3*i];        a1   = iatom[3*i+1];        a2   = iatom[3*i+2];        lenA = idef->iparams[type].constr.dA;        lenB = idef->iparams[type].constr.dB;        /* Skip the flexible constraints when not doing dynamics */        if (bDynamics || lenA!=0 || lenB!=0)        {            li->bllen0[con]  = lenA;            li->ddist[con]   = lenB - lenA;//.........这里部分代码省略.........

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

gmx_bool replica_exchange(FILE *fplog, const t_commrec *cr, struct gmx_repl_ex *re,                          t_state *state, gmx_enerdata_t *enerd,                          t_state *state_local, gmx_int64_t step, real time){    int i, j;    int replica_id = 0;    int exchange_partner;    int maxswap = 0;    /* Number of rounds of exchanges needed to deal with any multiple     * exchanges. */    /* Where each replica ends up after the exchange attempt(s). */    /* The order in which multiple exchanges will occur. */    gmx_bool bThisReplicaExchanged = FALSE;    if (MASTER(cr))    {        replica_id  = re->repl;        test_for_replica_exchange(fplog, cr->ms, re, enerd, det(state_local->box), step, time);        prepare_to_do_exchange(fplog, re->destinations, replica_id, re->nrepl, &maxswap,                               re->order, re->cyclic, re->incycle, &bThisReplicaExchanged);    }    /* Do intra-simulation broadcast so all processors belonging to     * each simulation know whether they need to participate in     * collecting the state. Otherwise, they might as well get on with     * the next thing to do. */    if (DOMAINDECOMP(cr))    {#ifdef GMX_MPI        MPI_Bcast(&bThisReplicaExchanged, sizeof(gmx_bool), MPI_BYTE, MASTERRANK(cr),                  cr->mpi_comm_mygroup);#endif    }    if (bThisReplicaExchanged)    {        /* Exchange the states */        /* Collect the global state on the master node */        if (DOMAINDECOMP(cr))        {            dd_collect_state(cr->dd, state_local, state);        }        else        {            copy_state_nonatomdata(state_local, state);        }        if (MASTER(cr))        {            /* There will be only one swap cycle with standard replica             * exchange, but there may be multiple swap cycles if we             * allow multiple swaps. */            for (j = 0; j < maxswap; j++)            {                exchange_partner = re->order[replica_id][j];                if (exchange_partner != replica_id)                {                    /* Exchange the global states between the master nodes */                    if (debug)                    {                        fprintf(debug, "Exchanging %d with %d/n", replica_id, exchange_partner);                    }                    exchange_state(cr->ms, exchange_partner, state);                }            }            /* For temperature-type replica exchange, we need to scale             * the velocities. */            if (re->type == ereTEMP || re->type == ereTL)            {                scale_velocities(state, sqrt(re->q[ereTEMP][replica_id]/re->q[ereTEMP][re->destinations[replica_id]]));            }        }        /* With domain decomposition the global state is distributed later */        if (!DOMAINDECOMP(cr))        {            /* Copy the global state to the local state data structure */            copy_state_nonatomdata(state, state_local);        }    }    return bThisReplicaExchanged;}

void mdAlgorithmsSetupAtomData(t_commrec         *cr,                               const t_inputrec  *ir,                               const gmx_mtop_t  *top_global,                               gmx_localtop_t    *top,                               t_forcerec        *fr,                               t_graph          **graph,                               t_mdatoms         *mdatoms,                               gmx_vsite_t       *vsite,                               gmx_shellfc_t     *shellfc){    bool  usingDomDec = DOMAINDECOMP(cr);    int   numAtomIndex, numHomeAtoms;    int  *atomIndex;    if (usingDomDec)    {        numAtomIndex = dd_natoms_mdatoms(cr->dd);        atomIndex    = cr->dd->gatindex;        numHomeAtoms = cr->dd->nat_home;    }    else    {        numAtomIndex = -1;        atomIndex    = NULL;        numHomeAtoms = top_global->natoms;    }    atoms2md(top_global, ir, numAtomIndex, atomIndex, numHomeAtoms, mdatoms);    if (usingDomDec)    {        dd_sort_local_top(cr->dd, mdatoms, top);    }    else    {        /* Currently gmx_generate_local_top allocates and returns a pointer.         * We should implement a more elegant solution.         */        gmx_localtop_t *tmpTop;        tmpTop = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);        *top   = *tmpTop;        sfree(tmpTop);    }    if (vsite)    {        if (usingDomDec)        {            /* The vsites were already assigned by the domdec topology code.             * We only need to do the thread division here.             */            split_vsites_over_threads(top->idef.il, top->idef.iparams,                                      mdatoms, FALSE, vsite);        }        else        {            set_vsite_top(vsite, top, mdatoms, cr);        }    }    if (!usingDomDec && ir->ePBC != epbcNONE && !fr->bMolPBC)    {        GMX_ASSERT(graph != NULL, "We use a graph with PBC (no periodic mols) and without DD");        *graph = mk_graph(NULL, &(top->idef), 0, top_global->natoms, FALSE, FALSE);    }    else if (graph != NULL)    {        *graph = NULL;    }    /* Note that with DD only flexible constraints, not shells, are supported     * and these don't require setup in make_local_shells().     */    if (!usingDomDec && shellfc)    {        make_local_shells(cr, mdatoms, shellfc);    }    setup_bonded_threading(fr, &top->idef);}

//.........这里部分代码省略.........            pme_lb->rbuf_vdw = ic->rlistlong - ic->rvdw;        }        else        {            pme_lb->rbuf_vdw = ic->rlist - ic->rvdw;        }    }    copy_mat(box, pme_lb->box_start);    if (ir->ePBC == epbcXY && ir->nwall == 2)    {        svmul(ir->wall_ewald_zfac, pme_lb->box_start[ZZ], pme_lb->box_start[ZZ]);    }    pme_lb->n = 1;    snew(pme_lb->setup, pme_lb->n);    pme_lb->rcut_vdw                 = ic->rvdw;    pme_lb->rcut_coulomb_start       = ir->rcoulomb;    pme_lb->nstcalclr_start          = ir->nstcalclr;    pme_lb->cur                      = 0;    pme_lb->setup[0].rcut_coulomb    = ic->rcoulomb;    pme_lb->setup[0].rlist           = ic->rlist;    pme_lb->setup[0].rlistlong       = ic->rlistlong;    pme_lb->setup[0].nstcalclr       = ir->nstcalclr;    pme_lb->setup[0].grid[XX]        = ir->nkx;    pme_lb->setup[0].grid[YY]        = ir->nky;    pme_lb->setup[0].grid[ZZ]        = ir->nkz;    pme_lb->setup[0].ewaldcoeff_q    = ic->ewaldcoeff_q;    pme_lb->setup[0].ewaldcoeff_lj   = ic->ewaldcoeff_lj;    pme_lb->setup[0].pmedata         = pmedata;    spm = 0;    for (d = 0; d < DIM; d++)    {        sp = norm(pme_lb->box_start[d])/pme_lb->setup[0].grid[d];        if (sp > spm)        {            spm = sp;        }    }    pme_lb->setup[0].spacing = spm;    if (ir->fourier_spacing > 0)    {        pme_lb->cut_spacing = ir->rcoulomb/ir->fourier_spacing;    }    else    {        pme_lb->cut_spacing = ir->rcoulomb/pme_lb->setup[0].spacing;    }    pme_lb->stage = 0;    pme_lb->fastest     = 0;    pme_lb->lower_limit = 0;    pme_lb->start       = 0;    pme_lb->end         = 0;    pme_lb->elimited    = epmelblimNO;    pme_lb->cycles_n = 0;    pme_lb->cycles_c = 0;    /* Tune with GPUs and/or separate PME ranks.     * When running only on a CPU without PME ranks, PME tuning will only help     * with small numbers of atoms in the cut-off sphere.     */    pme_lb->bActive  = (wallcycle_have_counter() && (bUseGPU ||                                                     pme_lb->bSepPMERanks));    /* With GPUs and no separate PME ranks we can't measure the PP/PME     * imbalance, so we start balancing right away.     * Otherwise we only start balancing after we observe imbalance.     */    pme_lb->bBalance = (pme_lb->bActive && (bUseGPU && !pme_lb->bSepPMERanks));    pme_lb->step_rel_stop = PMETunePeriod*ir->nstlist;    /* Delay DD load balancing when GPUs are used */    if (pme_lb->bActive && DOMAINDECOMP(cr) && cr->dd->nnodes > 1 && bUseGPU)    {        /* Lock DLB=auto to off (does nothing when DLB=yes/no.         * With GPUs and separate PME nodes, we want to first         * do PME tuning without DLB, since DLB might limit         * the cut-off, which never improves performance.         * We allow for DLB + PME tuning after a first round of tuning.         */        dd_dlb_lock(cr->dd);        if (dd_dlb_is_locked(cr->dd))        {            md_print_warn(cr, fp_log, "NOTE: DLB will not turn on during the first phase of PME tuning/n");        }    }    *pme_lb_p = pme_lb;    *bPrinting = pme_lb->bBalance;}