double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
t_trxstatus *status;
t_trxframe rerun_fr;
- gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS, bOurStep;
+ gmx_bool bDispCorr, bCharge, bRFExcl, bNotLastFrame, bStateChanged, bNS;
tensor force_vir, shake_vir, vir, pres;
int cg_tp, a_tp0, a_tp1, ngid, gid_tp, nener, e;
rvec *x_mol;
rvec mu_tot, x_init, dx, x_tp;
- int nnodes, frame, nsteps, step;
+ int nnodes, frame;
+ gmx_int64_t frame_step_prev, frame_step;
+ gmx_int64_t nsteps, stepblocksize = 0, step;
+ gmx_int64_t rnd_count_stride, rnd_count;
+ gmx_int64_t seed;
+ double rnd[4];
int i, start, end;
- gmx_rng_t tpi_rand;
FILE *fp_tpi = NULL;
char *ptr, *dump_pdb, **leg, str[STRLEN], str2[STRLEN];
double dbl, dump_ener;
a_tp1-a_tp0, bCharge ? "with" : "without");
fprintf(fplog, "\nWill insert %d times in each frame of %s\n",
- nsteps, opt2fn("-rerun", nfile, fnm));
+ (int)nsteps, opt2fn("-rerun", nfile, fnm));
}
if (!bCavity)
}
snew(sum_UgembU, nener);
- /* Initialize random generator */
- tpi_rand = gmx_rng_init(inputrec->ld_seed);
+ /* Copy the random seed set by the user */
+ seed = inputrec->ld_seed;
+ /* We use the frame step number as one random counter.
+ * The second counter use the insertion (step) count. But we
+ * need multiple random numbers per insertion. This number is
+ * not fixed, since we generate random locations in a sphere
+ * by putting locations in a cube and some of these fail.
+ * A count of 20 is already extremely unlikely, so 10000 is
+ * a safe margin for random numbers per insertion.
+ */
+ rnd_count_stride = 10000;
if (MASTER(cr))
{
nbin = 10;
snew(bin, nbin);
+ /* Avoid frame step numbers <= -1 */
+ frame_step_prev = -1;
+
bNotLastFrame = read_first_frame(oenv, &status, opt2fn("-rerun", nfile, fnm),
&rerun_fr, TRX_NEED_X);
frame = 0;
refvolshift = log(det(rerun_fr.box));
+ switch (inputrec->eI)
+ {
+ case eiTPI:
+ stepblocksize = inputrec->nstlist;
+ break;
+ case eiTPIC:
+ stepblocksize = 1;
+ break;
+ default:
+ gmx_fatal(FARGS, "Unknown integrator %s", ei_names[inputrec->eI]);
+ }
+
#ifdef GMX_SIMD_X86_SSE2_OR_HIGHER
/* Make sure we don't detect SSE overflow generated before this point */
gmx_mm_check_and_reset_overflow();
while (bNotLastFrame)
{
+ frame_step = rerun_fr.step;
+ if (frame_step <= frame_step_prev)
+ {
+ /* We don't have step number in the trajectory file,
+ * or we have constant or decreasing step numbers.
+ * Ensure we have increasing step numbers, since we use
+ * the step numbers as a counter for random numbers.
+ */
+ frame_step = frame_step_prev + 1;
+ }
+ frame_step_prev = frame_step;
+
lambda = rerun_fr.lambda;
t = rerun_fr.time;
bStateChanged = TRUE;
bNS = TRUE;
- for (step = 0; step < nsteps; step++)
+
+ step = cr->nodeid*stepblocksize;
+ while (step < nsteps)
{
- /* In parallel all nodes generate all random configurations.
- * In that way the result is identical to a single cpu tpi run.
+ /* Initialize the second counter for random numbers using
+ * the insertion step index. This ensures that we get
+ * the same random numbers independently of how many
+ * MPI ranks we use. Also for the same seed, we get
+ * the same initial random sequence for different nsteps.
*/
+ rnd_count = step*rnd_count_stride;
+
if (!bCavity)
{
/* Random insertion in the whole volume */
if (bNS)
{
/* Generate a random position in the box */
- x_init[XX] = gmx_rng_uniform_real(tpi_rand)*state->box[XX][XX];
- x_init[YY] = gmx_rng_uniform_real(tpi_rand)*state->box[YY][YY];
- x_init[ZZ] = gmx_rng_uniform_real(tpi_rand)*state->box[ZZ][ZZ];
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd);
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd+2);
+ for (d = 0; d < DIM; d++)
+ {
+ x_init[d] = rnd[d]*state->box[d][d];
+ }
}
if (inputrec->nstlist == 1)
{
/* Generate coordinates within |dx|=drmax of x_init */
do
{
- dx[XX] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
- dx[YY] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
- dx[ZZ] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd);
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd+2);
+ for (d = 0; d < DIM; d++)
+ {
+ dx[d] = (2*rnd[d] - 1)*drmax;
+ }
}
while (norm2(dx) > drmax*drmax);
rvec_add(x_init, dx, x_tp);
/* Generate coordinates within |dx|=drmax of x_init */
do
{
- dx[XX] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
- dx[YY] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
- dx[ZZ] = (2*gmx_rng_uniform_real(tpi_rand) - 1)*drmax;
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd);
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd+2);
+ for (d = 0; d < DIM; d++)
+ {
+ dx[d] = (2*rnd[d] - 1)*drmax;
+ }
}
while (norm2(dx) > drmax*drmax);
rvec_add(x_init, dx, x_tp);
copy_rvec(x_mol[i-a_tp0], state->x[i]);
}
/* Rotate the molecule randomly */
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd);
+ gmx_rng_cycle_2uniform(frame_step, rnd_count++, seed, RND_SEED_TPI, rnd+2);
rotate_conf(a_tp1-a_tp0, state->x+a_tp0, NULL,
- 2*M_PI*gmx_rng_uniform_real(tpi_rand),
- 2*M_PI*gmx_rng_uniform_real(tpi_rand),
- 2*M_PI*gmx_rng_uniform_real(tpi_rand));
+ 2*M_PI*rnd[0],
+ 2*M_PI*rnd[1],
+ 2*M_PI*rnd[2]);
/* Shift to the insertion location */
for (i = a_tp0; i < a_tp1; i++)
{
}
}
- /* Check if this insertion belongs to this node */
- bOurStep = TRUE;
- if (PAR(cr))
+ /* Clear some matrix variables */
+ clear_mat(force_vir);
+ clear_mat(shake_vir);
+ clear_mat(vir);
+ clear_mat(pres);
+
+ /* Set the charge group center of mass of the test particle */
+ copy_rvec(x_init, fr->cg_cm[top->cgs.nr-1]);
+
+ /* Calc energy (no forces) on new positions.
+ * Since we only need the intermolecular energy
+ * and the RF exclusion terms of the inserted molecule occur
+ * within a single charge group we can pass NULL for the graph.
+ * This also avoids shifts that would move charge groups
+ * out of the box.
+ *
+ * Some checks above ensure than we can not have
+ * twin-range interactions together with nstlist > 1,
+ * therefore we do not need to remember the LR energies.
+ */
+ /* Make do_force do a single node force calculation */
+ cr->nnodes = 1;
+ do_force(fplog, cr, inputrec,
+ step, nrnb, wcycle, top, &top_global->groups,
+ state->box, state->x, &state->hist,
+ f, force_vir, mdatoms, enerd, fcd,
+ state->lambda,
+ NULL, fr, NULL, mu_tot, t, NULL, NULL, FALSE,
+ GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
+ (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS | GMX_FORCE_DO_LR : 0) |
+ (bStateChanged ? GMX_FORCE_STATECHANGED : 0));
+ cr->nnodes = nnodes;
+ bStateChanged = FALSE;
+ bNS = FALSE;
+
+ /* Calculate long range corrections to pressure and energy */
+ calc_dispcorr(fplog, inputrec, fr, step, top_global->natoms, state->box,
+ lambda, pres, vir, &prescorr, &enercorr, &dvdlcorr);
+ /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
+ enerd->term[F_DISPCORR] = enercorr;
+ enerd->term[F_EPOT] += enercorr;
+ enerd->term[F_PRES] += prescorr;
+ enerd->term[F_DVDL_VDW] += dvdlcorr;
+
+ epot = enerd->term[F_EPOT];
+ bEnergyOutOfBounds = FALSE;
+#ifdef GMX_SIMD_X86_SSE2_OR_HIGHER
+ /* With SSE the energy can overflow, check for this */
+ if (gmx_mm_check_and_reset_overflow())
{
- switch (inputrec->eI)
+ if (debug)
{
- case eiTPI:
- bOurStep = ((step / inputrec->nstlist) % nnodes == cr->nodeid);
- break;
- case eiTPIC:
- bOurStep = (step % nnodes == cr->nodeid);
- break;
- default:
- gmx_fatal(FARGS, "Unknown integrator %s", ei_names[inputrec->eI]);
+ fprintf(debug, "Found an SSE overflow, assuming the energy is out of bounds\n");
}
+ bEnergyOutOfBounds = TRUE;
}
- if (bOurStep)
- {
- /* Clear some matrix variables */
- clear_mat(force_vir);
- clear_mat(shake_vir);
- clear_mat(vir);
- clear_mat(pres);
-
- /* Set the charge group center of mass of the test particle */
- copy_rvec(x_init, fr->cg_cm[top->cgs.nr-1]);
-
- /* Calc energy (no forces) on new positions.
- * Since we only need the intermolecular energy
- * and the RF exclusion terms of the inserted molecule occur
- * within a single charge group we can pass NULL for the graph.
- * This also avoids shifts that would move charge groups
- * out of the box.
- *
- * Some checks above ensure than we can not have
- * twin-range interactions together with nstlist > 1,
- * therefore we do not need to remember the LR energies.
- */
- /* Make do_force do a single node force calculation */
- cr->nnodes = 1;
- do_force(fplog, cr, inputrec,
- step, nrnb, wcycle, top, &top_global->groups,
- state->box, state->x, &state->hist,
- f, force_vir, mdatoms, enerd, fcd,
- state->lambda,
- NULL, fr, NULL, mu_tot, t, NULL, NULL, FALSE,
- GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
- (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS | GMX_FORCE_DO_LR : 0) |
- (bStateChanged ? GMX_FORCE_STATECHANGED : 0));
- cr->nnodes = nnodes;
- bStateChanged = FALSE;
- bNS = FALSE;
-
- /* Calculate long range corrections to pressure and energy */
- calc_dispcorr(fplog, inputrec, fr, step, top_global->natoms, state->box,
- lambda, pres, vir, &prescorr, &enercorr, &dvdlcorr);
- /* figure out how to rearrange the next 4 lines MRS 8/4/2009 */
- enerd->term[F_DISPCORR] = enercorr;
- enerd->term[F_EPOT] += enercorr;
- enerd->term[F_PRES] += prescorr;
- enerd->term[F_DVDL_VDW] += dvdlcorr;
-
- epot = enerd->term[F_EPOT];
- bEnergyOutOfBounds = FALSE;
-#ifdef GMX_SIMD_X86_SSE2_OR_HIGHER
- /* With SSE the energy can overflow, check for this */
- if (gmx_mm_check_and_reset_overflow())
- {
- if (debug)
- {
- fprintf(debug, "Found an SSE overflow, assuming the energy is out of bounds\n");
- }
- bEnergyOutOfBounds = TRUE;
- }
#endif
- /* If the compiler doesn't optimize this check away
- * we catch the NAN energies.
- * The epot>GMX_REAL_MAX check catches inf values,
- * which should nicely result in embU=0 through the exp below,
- * but it does not hurt to check anyhow.
- */
- /* Non-bonded Interaction usually diverge at r=0.
- * With tabulated interaction functions the first few entries
- * should be capped in a consistent fashion between
- * repulsion, dispersion and Coulomb to avoid accidental
- * negative values in the total energy.
- * The table generation code in tables.c does this.
- * With user tbales the user should take care of this.
- */
- if (epot != epot || epot > GMX_REAL_MAX)
+ /* If the compiler doesn't optimize this check away
+ * we catch the NAN energies.
+ * The epot>GMX_REAL_MAX check catches inf values,
+ * which should nicely result in embU=0 through the exp below,
+ * but it does not hurt to check anyhow.
+ */
+ /* Non-bonded Interaction usually diverge at r=0.
+ * With tabulated interaction functions the first few entries
+ * should be capped in a consistent fashion between
+ * repulsion, dispersion and Coulomb to avoid accidental
+ * negative values in the total energy.
+ * The table generation code in tables.c does this.
+ * With user tbales the user should take care of this.
+ */
+ if (epot != epot || epot > GMX_REAL_MAX)
+ {
+ bEnergyOutOfBounds = TRUE;
+ }
+ if (bEnergyOutOfBounds)
+ {
+ if (debug)
{
- bEnergyOutOfBounds = TRUE;
+ fprintf(debug, "\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t, (int)step, epot);
}
- if (bEnergyOutOfBounds)
+ embU = 0;
+ }
+ else
+ {
+ embU = exp(-beta*epot);
+ sum_embU += embU;
+ /* Determine the weighted energy contributions of each energy group */
+ e = 0;
+ sum_UgembU[e++] += epot*embU;
+ if (fr->bBHAM)
{
- if (debug)
+ for (i = 0; i < ngid; i++)
{
- fprintf(debug, "\n time %.3f, step %d: non-finite energy %f, using exp(-bU)=0\n", t, step, epot);
+ sum_UgembU[e++] +=
+ (enerd->grpp.ener[egBHAMSR][GID(i, gid_tp, ngid)] +
+ enerd->grpp.ener[egBHAMLR][GID(i, gid_tp, ngid)])*embU;
}
- embU = 0;
}
else
{
- embU = exp(-beta*epot);
- sum_embU += embU;
- /* Determine the weighted energy contributions of each energy group */
- e = 0;
- sum_UgembU[e++] += epot*embU;
- if (fr->bBHAM)
- {
- for (i = 0; i < ngid; i++)
- {
- sum_UgembU[e++] +=
- (enerd->grpp.ener[egBHAMSR][GID(i, gid_tp, ngid)] +
- enerd->grpp.ener[egBHAMLR][GID(i, gid_tp, ngid)])*embU;
- }
- }
- else
- {
- for (i = 0; i < ngid; i++)
- {
- sum_UgembU[e++] +=
- (enerd->grpp.ener[egLJSR][GID(i, gid_tp, ngid)] +
- enerd->grpp.ener[egLJLR][GID(i, gid_tp, ngid)])*embU;
- }
- }
- if (bDispCorr)
+ for (i = 0; i < ngid; i++)
{
- sum_UgembU[e++] += enerd->term[F_DISPCORR]*embU;
- }
- if (bCharge)
- {
- for (i = 0; i < ngid; i++)
- {
- sum_UgembU[e++] +=
- (enerd->grpp.ener[egCOULSR][GID(i, gid_tp, ngid)] +
- enerd->grpp.ener[egCOULLR][GID(i, gid_tp, ngid)])*embU;
- }
- if (bRFExcl)
- {
- sum_UgembU[e++] += enerd->term[F_RF_EXCL]*embU;
- }
- if (EEL_FULL(fr->eeltype))
- {
- sum_UgembU[e++] += enerd->term[F_COUL_RECIP]*embU;
- }
+ sum_UgembU[e++] +=
+ (enerd->grpp.ener[egLJSR][GID(i, gid_tp, ngid)] +
+ enerd->grpp.ener[egLJLR][GID(i, gid_tp, ngid)])*embU;
}
}
-
- if (embU == 0 || beta*epot > bU_bin_limit)
+ if (bDispCorr)
{
- bin[0]++;
+ sum_UgembU[e++] += enerd->term[F_DISPCORR]*embU;
}
- else
+ if (bCharge)
{
- i = (int)((bU_logV_bin_limit
- - (beta*epot - logV + refvolshift))*invbinw
- + 0.5);
- if (i < 0)
+ for (i = 0; i < ngid; i++)
+ {
+ sum_UgembU[e++] +=
+ (enerd->grpp.ener[egCOULSR][GID(i, gid_tp, ngid)] +
+ enerd->grpp.ener[egCOULLR][GID(i, gid_tp, ngid)])*embU;
+ }
+ if (bRFExcl)
{
- i = 0;
+ sum_UgembU[e++] += enerd->term[F_RF_EXCL]*embU;
}
- if (i >= nbin)
+ if (EEL_FULL(fr->eeltype))
{
- realloc_bins(&bin, &nbin, i+10);
+ sum_UgembU[e++] += enerd->term[F_COUL_RECIP]*embU;
}
- bin[i]++;
}
+ }
- if (debug)
+ if (embU == 0 || beta*epot > bU_bin_limit)
+ {
+ bin[0]++;
+ }
+ else
+ {
+ i = (int)((bU_logV_bin_limit
+ - (beta*epot - logV + refvolshift))*invbinw
+ + 0.5);
+ if (i < 0)
{
- fprintf(debug, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
- step, epot, x_tp[XX], x_tp[YY], x_tp[ZZ]);
+ i = 0;
}
-
- if (dump_pdb && epot <= dump_ener)
+ if (i >= nbin)
{
- sprintf(str, "t%g_step%d.pdb", t, step);
- sprintf(str2, "t: %f step %d ener: %f", t, step, epot);
- write_sto_conf_mtop(str, str2, top_global, state->x, state->v,
- inputrec->ePBC, state->box);
+ realloc_bins(&bin, &nbin, i+10);
}
+ bin[i]++;
+ }
+
+ if (debug)
+ {
+ fprintf(debug, "TPI %7d %12.5e %12.5f %12.5f %12.5f\n",
+ (int)step, epot, x_tp[XX], x_tp[YY], x_tp[ZZ]);
+ }
+
+ if (dump_pdb && epot <= dump_ener)
+ {
+ sprintf(str, "t%g_step%d.pdb", t, (int)step);
+ sprintf(str2, "t: %f step %d ener: %f", t, (int)step, epot);
+ write_sto_conf_mtop(str, str2, top_global, state->x, state->v,
+ inputrec->ePBC, state->box);
+ }
+
+ step++;
+ if ((step/stepblocksize) % cr->nnodes != cr->nodeid)
+ {
+ /* Skip all steps assigned to the other MPI ranks */
+ step += (cr->nnodes - 1)*stepblocksize;
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff