X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fmdlib%2Fcoupling.c;h=b156c61fad22009667945749571563f0af1b46bd;hb=e336499626008e8ff06b1648c1b7afbab38c949b;hp=3d1190e401e6a87a33f0edd9add447de043c59a7;hpb=16200c9bc8d86775ddc0f8823f663a0eb7944b50;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/mdlib/coupling.c b/src/gromacs/mdlib/coupling.c index 3d1190e401..b156c61fad 100644 --- a/src/gromacs/mdlib/coupling.c +++ b/src/gromacs/mdlib/coupling.c @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013, by the GROMACS development team, led by + * Copyright (c) 2013,2014, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -34,24 +34,23 @@ * To help us fund GROMACS development, we humbly ask that you cite * the research papers on the package. Check out http://www.gromacs.org. */ -#ifdef HAVE_CONFIG_H -#include -#endif +#include "gmxpre.h" + #include -#include "typedefs.h" -#include "smalloc.h" -#include "update.h" -#include "vec.h" -#include "macros.h" -#include "physics.h" -#include "names.h" -#include "gmx_fatal.h" -#include "txtdump.h" -#include "nrnb.h" -#include "gmx_random.h" -#include "update.h" -#include "mdrun.h" +#include "gromacs/legacyheaders/macros.h" +#include "gromacs/legacyheaders/mdrun.h" +#include "gromacs/legacyheaders/names.h" +#include "gromacs/legacyheaders/nrnb.h" +#include "gromacs/legacyheaders/txtdump.h" +#include "gromacs/legacyheaders/typedefs.h" +#include "gromacs/legacyheaders/update.h" +#include "gromacs/legacyheaders/types/commrec.h" +#include "gromacs/math/units.h" +#include "gromacs/math/vec.h" +#include "gromacs/random/random.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/smalloc.h" #define NTROTTERPARTS 3 @@ -722,203 +721,50 @@ void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt) } } -static int poisson_variate(real lambda, gmx_rng_t rng) +void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step, + const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac) { + const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL); + int i; + int gc = 0; - real L; - int k = 0; - real p = 1.0; - - L = exp(-lambda); - - do - { - k = k+1; - p *= gmx_rng_uniform_real(rng); - } - while (p > L); - - return k-1; -} - -void andersen_tcoupl(t_inputrec *ir, t_mdatoms *md, t_state *state, gmx_rng_t rng, real rate, t_idef *idef, int nblocks, int *sblock, gmx_bool *randatom, int *randatom_list, gmx_bool *randomize, real *boltzfac) -{ - t_grpopts *opts; - int i, j, k, d, len, n, ngtc, gc = 0; - int nshake, nsettle, nrandom, nrand_group; - real boltz, scal, reft, prand; - t_iatom *iatoms; - - /* convenience variables */ - opts = &ir->opts; - ngtc = opts->ngtc; + /* randomize the velocities of the selected particles */ - /* idef is only passed in if it's chance-per-particle andersen, so - it essentially serves as a boolean to determine which type of - andersen is being used */ - if (idef) + for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */ { + int ng = gatindex ? gatindex[i] : i; + gmx_bool bRandomize; - /* randomly atoms to randomize. However, all constraint - groups have to have either all of the atoms or none of the - atoms randomize. - - Algorithm: - 1. Select whether or not to randomize each atom to get the correct probability. - 2. Cycle through the constraint groups. - 2a. for each constraint group, determine the fraction f of that constraint group that are - chosen to be randomized. - 2b. all atoms in the constraint group are randomized with probability f. - */ - - nrandom = 0; - if ((rate < 0.05) && (md->homenr > 50)) - { - /* if the rate is relatively high, use a standard method, if low rate, - * use poisson */ - /* poisson distributions approxmation, more efficient for - * low rates, fewer random numbers required */ - nrandom = poisson_variate(md->homenr*rate, rng); /* how many do we randomize? Use Poisson. */ - /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible. - worst thing that happens, it lowers the true rate an negligible amount */ - for (i = 0; i < nrandom; i++) - { - randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE; - } - } - else + if (md->cTC) { - for (i = 0; i < md->homenr; i++) - { - if (gmx_rng_uniform_real(rng) < rate) - { - randatom[i] = TRUE; - nrandom++; - } - } + gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */ } - - /* instead of looping over the constraint groups, if we had a - list of which atoms were in which constraint groups, we - could then loop over only the groups that are randomized - now. But that is not available now. Create later after - determining whether there actually is any slowing. */ - - /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */ - - nsettle = idef->il[F_SETTLE].nr/2; - for (i = 0; i < nsettle; i++) + if (randomize[gc]) { - iatoms = idef->il[F_SETTLE].iatoms; - nrand_group = 0; - for (k = 0; k < 3; k++) /* settles are always 3 atoms, hardcoded */ + if (ir->etc == etcANDERSEN) { - if (randatom[iatoms[2*i+1]+k]) - { - nrand_group++; /* count the number of atoms to be shaken in the settles group */ - randatom[iatoms[2*i+1]+k] = FALSE; - nrandom--; - } + bRandomize = TRUE; } - if (nrand_group > 0) + else { - prand = (nrand_group)/3.0; /* use this fraction to compute the probability the - whole group is randomized */ - if (gmx_rng_uniform_real(rng) < prand) - { - for (k = 0; k < 3; k++) - { - randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */ - } - nrandom += 3; - } - } - } + double uniform[2]; - /* now loop through the shake groups */ - nshake = nblocks; - for (i = 0; i < nshake; i++) - { - iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]); - len = sblock[i+1]-sblock[i]; - nrand_group = 0; - for (k = 0; k < len; k++) - { - if (k%3 != 0) - { /* only 2/3 of the sblock items are atoms, the others are labels */ - if (randatom[iatoms[k]]) - { - nrand_group++; - randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than - one group in the shake block */ - nrandom--; - } - } + gmx_rng_cycle_2uniform(step*2, ng, ir->andersen_seed, RND_SEED_ANDERSEN, uniform); + bRandomize = (uniform[0] < rate); } - if (nrand_group > 0) + if (bRandomize) { - prand = (nrand_group)/(1.0*(2*len/3)); - if (gmx_rng_uniform_real(rng) < prand) - { - for (k = 0; k < len; k++) - { - if (k%3 != 0) - { /* only 2/3 of the sblock items are atoms, the others are labels */ - randatom[iatoms[k]] = TRUE; /* randomize all of them */ - nrandom++; - } - } - } - } - } - if (nrandom > 0) - { - n = 0; - for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */ - { - if (randatom[i]) + real scal, gauss[3]; + int d; + + scal = sqrt(boltzfac[gc]*md->invmass[i]); + gmx_rng_cycle_3gaussian_table(step*2+1, ng, ir->andersen_seed, RND_SEED_ANDERSEN, gauss); + for (d = 0; d < DIM; d++) { - randatom_list[n] = i; - n++; + state->v[i][d] = scal*gauss[d]; } } - nrandom = n; /* there are some values of nrandom for which - this algorithm won't work; for example all - water molecules and nrandom =/= 3. Better to - recount and use this number (which we - calculate anyway: it will not affect - the average number of atoms accepted. - */ - } - } - else - { - /* if it's andersen-massive, then randomize all the atoms */ - nrandom = md->homenr; - for (i = 0; i < nrandom; i++) - { - randatom_list[i] = i; - } - } - - /* randomize the velocities of the selected particles */ - - for (i = 0; i < nrandom; i++) /* now loop over the list of atoms */ - { - n = randatom_list[i]; - if (md->cTC) - { - gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */ - } - if (randomize[gc]) - { - scal = sqrt(boltzfac[gc]*md->invmass[n]); - for (d = 0; d < DIM; d++) - { - state->v[n][d] = scal*gmx_rng_gaussian_table(rng); - } } - randatom[n] = FALSE; /* unmark this atom for randomization */ } } @@ -1059,7 +905,7 @@ void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind, /* but do we actually need the total? */ /* modify the velocities as well */ - for (n = md->start; n < md->start+md->homenr; n++) + for (n = 0; n < md->homenr; n++) { if (md->cTC) /* does this conditional need to be here? is this always true?*/ { @@ -1390,7 +1236,8 @@ int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool b real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ) { - int i, j, nd, ndj, bmass, qmass, ngtcall; + int i, j, bmass, qmass, ngtcall; + real nd, ndj; real ener_npt, reft, eta, kT, tau; double *ivxi, *ixi; double *iQinv; @@ -1473,7 +1320,7 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ) reft = max(ir->opts.ref_t[i], 0); kT = BOLTZ * reft; - if (nd > 0) + if (nd > 0.0) { if (IR_NVT_TROTTER(ir)) { @@ -1490,7 +1337,7 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ) } else { - ndj = 1; + ndj = 1.0; } ener_npt += ndj*ixi[j]*kT; } @@ -1507,86 +1354,109 @@ real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ) return ener_npt; } -static real vrescale_gamdev(int ia, gmx_rng_t rng) +static real vrescale_gamdev(real ia, + gmx_int64_t step, gmx_int64_t *count, + gmx_int64_t seed1, gmx_int64_t seed2) /* Gamma distribution, adapted from numerical recipes */ { - int j; - real am, e, s, v1, v2, x, y; + real am, e, s, v1, v2, x, y; + double rnd[2]; - if (ia < 6) - { - do - { - x = 1.0; - for (j = 1; j <= ia; j++) - { - x *= gmx_rng_uniform_real(rng); - } - } - while (x == 0); - x = -log(x); - } - else + assert(ia > 1); + + do { do { do { - do - { - v1 = gmx_rng_uniform_real(rng); - v2 = 2.0*gmx_rng_uniform_real(rng)-1.0; - } - while (v1*v1 + v2*v2 > 1.0 || - v1*v1*GMX_REAL_MAX < 3.0*ia); - /* The last check above ensures that both x (3.0 > 2.0 in s) - * and the pre-factor for e do not go out of range. - */ - y = v2/v1; - am = ia - 1; - s = sqrt(2.0*am + 1.0); - x = s*y + am; + gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd); + v1 = rnd[0]; + v2 = 2.0*rnd[1] - 1.0; } - while (x <= 0.0); - e = (1.0 + y*y)*exp(am*log(x/am) - s*y); + while (v1*v1 + v2*v2 > 1.0 || + v1*v1*GMX_REAL_MAX < 3.0*ia); + /* The last check above ensures that both x (3.0 > 2.0 in s) + * and the pre-factor for e do not go out of range. + */ + y = v2/v1; + am = ia - 1; + s = sqrt(2.0*am + 1.0); + x = s*y + am; } - while (gmx_rng_uniform_real(rng) > e); + while (x <= 0.0); + + e = (1.0 + y*y)*exp(am*log(x/am) - s*y); + + gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd); } + while (rnd[0] > e); return x; } -static real vrescale_sumnoises(int nn, gmx_rng_t rng) +static real gaussian_count(gmx_int64_t step, gmx_int64_t *count, + gmx_int64_t seed1, gmx_int64_t seed2) +{ + double rnd[2], x, y, r; + + do + { + gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd); + x = 2.0*rnd[0] - 1.0; + y = 2.0*rnd[1] - 1.0; + r = x*x + y*y; + } + while (r > 1.0 || r == 0.0); + + r = sqrt(-2.0*log(r)/r); + + return x*r; +} + +static real vrescale_sumnoises(real nn, + gmx_int64_t step, gmx_int64_t *count, + gmx_int64_t seed1, gmx_int64_t seed2) { /* - * Returns the sum of n independent gaussian noises squared + * Returns the sum of nn independent gaussian noises squared * (i.e. equivalent to summing the square of the return values - * of nn calls to gmx_rng_gaussian_real).xs + * of nn calls to gmx_rng_gaussian_real). */ - real rr; + const real ndeg_tol = 0.0001; + real r; - if (nn == 0) - { - return 0.0; - } - else if (nn == 1) + if (nn < 2 + ndeg_tol) { - rr = gmx_rng_gaussian_real(rng); - return rr*rr; - } - else if (nn % 2 == 0) - { - return 2.0*vrescale_gamdev(nn/2, rng); + int nn_int, i; + real gauss; + + nn_int = (int)(nn + 0.5); + + if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol) + { + gmx_fatal(FARGS, "The v-rescale thermostat was called with a group with #DOF=%f, but for #DOF<3 only integer #DOF are supported", nn + 1); + } + + r = 0; + for (i = 0; i < nn_int; i++) + { + gauss = gaussian_count(step, count, seed1, seed2); + + r += gauss*gauss; + } } else { - rr = gmx_rng_gaussian_real(rng); - return 2.0*vrescale_gamdev((nn-1)/2, rng) + rr*rr; + /* Use a gamma distribution for any real nn > 2 */ + r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2); } + + return r; } -static real vrescale_resamplekin(real kk, real sigma, int ndeg, real taut, - gmx_rng_t rng) +static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, + gmx_int64_t step, gmx_int64_t seed) { /* * Generates a new value for the kinetic energy, @@ -1596,7 +1466,11 @@ static real vrescale_resamplekin(real kk, real sigma, int ndeg, real taut, * ndeg: number of degrees of freedom of the atoms to be thermalized * taut: relaxation time of the thermostat, in units of 'how often this routine is called' */ - real factor, rr; + /* rnd_count tracks the step-local state for the cycle random + * number generator. + */ + gmx_int64_t rnd_count = 0; + real factor, rr, ekin_new; if (taut > 0.1) { @@ -1606,15 +1480,20 @@ static real vrescale_resamplekin(real kk, real sigma, int ndeg, real taut, { factor = 0.0; } - rr = gmx_rng_gaussian_real(rng); - return + + rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE); + + ekin_new = kk + - (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, rng) + rr*rr)/ndeg - kk) + + (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) + 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor); + + return ekin_new; } -void vrescale_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt, - double therm_integral[], gmx_rng_t rng) +void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step, + gmx_ekindata_t *ekind, real dt, + double therm_integral[]) { t_grpopts *opts; int i; @@ -1639,7 +1518,8 @@ void vrescale_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt, Ek_ref = Ek_ref1*opts->nrdf[i]; Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], - opts->tau_t[i]/dt, rng); + opts->tau_t[i]/dt, + step, ir->ld_seed); /* Analytically Ek_new>=0, but we check for rounding errors */ if (Ek_new <= 0)