From: Berk Hess Date: Sun, 1 Mar 2015 20:38:07 +0000 (+0100) Subject: Allow pull groups of 1 atom with mass 0 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=d9192a3e93e814902e50ac15c1552d8bb3c757cb;p=alexxy%2Fgromacs.git Allow pull groups of 1 atom with mass 0 Unless constraint pulling is used, a pull group consisting of 1 atom can have mass 0, since the mass of the COM is irrelevant. This is useful for pulling on a virtual site. Also reorganized conditional and loop orders in the COM calculation code and made many variables local to the scope they are used in. Change-Id: Ic2950d2db19df8673c0d2e5090ce6e121e45bf2d --- diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c index fca9f171c9..539d015e19 100644 --- a/src/gromacs/pulling/pull.c +++ b/src/gromacs/pulling/pull.c @@ -302,18 +302,32 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, inv_wm = pgrp->mwscale; - for (i = 0; i < pgrp->nat_loc; i++) + if (pgrp->nat == 1 && pgrp->nat_loc == 1) { - ii = pgrp->ind_loc[i]; - wmass = md->massT[ii]; - if (pgrp->weight_loc) + /* Only one atom and our rank has this atom: we can skip + * the mass weighting, which means that this code also works + * for mass=0, e.g. with a virtual site. + */ + for (m = 0; m < DIM; m++) { - wmass *= pgrp->weight_loc[i]; + f[pgrp->ind_loc[0]][m] += sign*f_pull[m]; } - - for (m = 0; m < DIM; m++) + } + else + { + for (i = 0; i < pgrp->nat_loc; i++) { - f[ii][m] += sign * wmass * f_pull[m] * inv_wm; + ii = pgrp->ind_loc[i]; + wmass = md->massT[ii]; + if (pgrp->weight_loc) + { + wmass *= pgrp->weight_loc[i]; + } + + for (m = 0; m < DIM; m++) + { + f[ii][m] += sign * wmass * f_pull[m] * inv_wm; + } } } } @@ -1219,8 +1233,19 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, if (wmass == 0) { - gmx_fatal(FARGS, "The total%s mass of pull group %d is zero", - pg->weight ? " weighted" : "", g); + /* We can have single atom groups with zero mass with potential pulling + * without cosine weighting. + */ + if (pg->nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS) + { + /* With one atom the mass doesn't matter */ + wwmass = 1; + } + else + { + gmx_fatal(FARGS, "The total%s mass of pull group %d is zero", + pg->weight ? " weighted" : "", g); + } } if (fplog) { diff --git a/src/gromacs/pulling/pullutil.c b/src/gromacs/pulling/pullutil.c index 5f231846dd..135edf649c 100644 --- a/src/gromacs/pulling/pullutil.c +++ b/src/gromacs/pulling/pullutil.c @@ -36,6 +36,7 @@ */ #include "gmxpre.h" +#include #include #include "gromacs/fileio/confio.h" @@ -324,13 +325,8 @@ 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; - dvec com, comp; - double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw; - rvec x_pbc = {0, 0, 0}, dx; - t_pull_group *pgrp; + int g; + real twopi_box = 0; if (pull->rbuf == NULL) { @@ -360,6 +356,10 @@ void pull_calc_coms(t_commrec *cr, 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) @@ -372,39 +372,51 @@ void pull_calc_coms(t_commrec *cr, for (g = 0; g < pull->ngroup; g++) { + t_pull_group *pgrp; + 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 (pgrp->bCalcCOM) { - 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++) + if (pgrp->epgrppbc != epgrppbcCOS) { - ii = pgrp->ind_loc[i]; - mass = md->massT[ii]; - 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(pull->rbuf[g], x_pbc); + } + + for (i = 0; i < pgrp->nat_loc; i++) { - if (pgrp->weight_loc) + int ii, m; + real mass, wm; + + ii = pgrp->ind_loc[i]; + mass = md->massT[ii]; + if (pgrp->weight_loc == NULL) { - w = pgrp->weight_loc[i]; + wm = mass; + wmass += wm; + } + else + { + real w; + + w = pgrp->weight_loc[i]; + wm = w*mass; + wmass += wm; + wwmass += wm*w; } - wm = w*mass; - wmass += wm; - wwmass += wm*w; if (pgrp->epgrppbc == epgrppbcNONE) { /* Plain COM: sum the coordinates */ @@ -422,6 +434,8 @@ void pull_calc_coms(t_commrec *cr, } else { + rvec dx; + /* Sum the difference with the reference atom */ pbc_dx(pbc, x[ii], x_pbc, dx); for (m = 0; m < DIM; m++) @@ -441,8 +455,60 @@ void pull_calc_coms(t_commrec *cr, } } } - else + + /* We do this check after the loop above to avoid more nesting. + * If we have a single-atom group the mass is irrelevant, so + * we can remove the mass factor to avoid division by zero. + * Note that with constraint pulling the mass does matter, but + * in that case a check group mass != 0 has been done before. + */ + if (pgrp->nat == 1 && pgrp->nat_loc == 1 && wmass == 0) { + int m; + + /* Copy the single atom coordinate */ + for (m = 0; m < DIM; m++) + { + com[m] = x[pgrp->ind_loc[0]][m]; + } + /* Set all mass factors to 1 to get the correct COM */ + wmass = 1; + wwmass = 1; + } + + if (pgrp->weight_loc == NULL) + { + wwmass = wmass; + } + + /* Copy local sums to a buffer for global summing */ + copy_dvec(com, pull->dbuf[g*3]); + copy_dvec(comp, pull->dbuf[g*3+1]); + pull->dbuf[g*3+2][0] = wmass; + pull->dbuf[g*3+2][1] = wwmass; + pull->dbuf[g*3+2][2] = 0; + } + else + { + /* Cosine weighting geometry */ + double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw; + int i; + + cm = 0; + sm = 0; + cmp = 0; + smp = 0; + ccm = 0; + csm = 0; + ssm = 0; + + for (i = 0; i < pgrp->nat_loc; i++) + { + int ii; + real mass; + + ii = pgrp->ind_loc[i]; + mass = md->massT[ii]; /* Determine cos and sin sums */ csw = cos(x[ii][pull->cosdim]*twopi_box); snw = sin(x[ii][pull->cosdim]*twopi_box); @@ -460,21 +526,8 @@ void pull_calc_coms(t_commrec *cr, smp += snw*mass; } } - } - } - /* Copy local sums to a buffer for global summing */ - switch (pgrp->epgrppbc) - { - case epgrppbcNONE: - case epgrppbcREFAT: - copy_dvec(com, pull->dbuf[g*3]); - copy_dvec(comp, pull->dbuf[g*3+1]); - pull->dbuf[g*3+2][0] = wmass; - pull->dbuf[g*3+2][1] = wwmass; - pull->dbuf[g*3+2][2] = 0; - break; - case epgrppbcCOS: + /* Copy local sums to a buffer for global summing */ pull->dbuf[g*3 ][0] = cm; pull->dbuf[g*3 ][1] = sm; pull->dbuf[g*3 ][2] = 0; @@ -484,7 +537,7 @@ void pull_calc_coms(t_commrec *cr, pull->dbuf[g*3+2][0] = cmp; pull->dbuf[g*3+2][1] = smp; pull->dbuf[g*3+2][2] = 0; - break; + } } } @@ -496,11 +549,16 @@ void pull_calc_coms(t_commrec *cr, for (g = 0; g < pull->ngroup; g++) { + t_pull_group *pgrp; + pgrp = &pull->group[g]; if (pgrp->nat > 0 && pgrp->bCalcCOM) { if (pgrp->epgrppbc != epgrppbcCOS) { + double wmass, wwmass; + int m; + /* Determine the inverse mass */ wmass = pull->dbuf[g*3+2][0]; wwmass = pull->dbuf[g*3+2][1]; @@ -531,6 +589,10 @@ void pull_calc_coms(t_commrec *cr, } else { + /* Cosine weighting geometry */ + double csw, snw, wmass, wwmass; + int i, ii; + /* Determine the optimal location of the cosine weight */ csw = pull->dbuf[g*3][0]; snw = pull->dbuf[g*3][1]; @@ -562,8 +624,8 @@ void pull_calc_coms(t_commrec *cr, } if (debug) { - fprintf(debug, "Pull group %d wmass %f wwmass %f invtm %f\n", - g, wmass, wwmass, pgrp->invtm); + fprintf(debug, "Pull group %d wmass %f invtm %f\n", + g, 1.0/pgrp->mwscale, pgrp->invtm); } } }