\item {\tt GMX_DLB_BASED_ON_FLOPS}: do domain-decomposition dynamic load balancing based on flop count rather than
measured time elapsed (default 0, meaning off).
This makes the load balancing reproducible, which can be useful for debugging purposes.
- A value of 1 uses the flops; a value > 1 adds (value - 1)*5\% of noise to the flops to increase the imbalance and the scaling.
+ A value of 1 uses the flops; a value $>1$ adds $(\mbox{value} - 1)\times5\%$ of noise to the flops to increase the imbalance and the scaling.
\item {\tt GMX_DLB_MAX_BOX_SCALING}: maximum percentage box scaling permitted per domain-decomposition
load-balancing step (default 10)
\item {\tt GMX_DD_RECORD_LOAD}: record DD load statistics for reporting at end of the run (default 1, meaning on)
{
corner[ZZ] = zones->size[z].x1[ZZ];
}
+ if (dd->ndim == 1 && dd->dim[0] < ZZ && ZZ < dd->npbcdim &&
+ box[ZZ][1 - dd->dim[0]] != 0)
+ {
+ /* With 1D domain decomposition the cg's are not in
+ * the triclinic box, but triclinic x-y and rectangular y/x-z.
+ * Shift the corner of the z-vector back to along the box
+ * vector of dimension d, so it will later end up at 0 along d.
+ * This can affect the location of this corner along dd->dim[0]
+ * through the matrix operation below if box[d][dd->dim[0]]!=0.
+ */
+ int d = 1 - dd->dim[0];
+
+ corner[d] -= corner[ZZ]*box[ZZ][d]/box[ZZ][ZZ];
+ }
/* Apply the triclinic couplings */
assert(ddbox->npbcdim <= DIM);
for (i = YY; i < ddbox->npbcdim; i++)
{
for (j = XX; j < i; j++)
{
- /* With 1D domain decomposition the cg's are not in
- * a triclinic box, but triclinic x-y and rectangular y/x-z.
- * So we should ignore the coupling for the non
- * domain-decomposed dimension of the pair x and y.
- */
- if (!(dd->ndim == 1 && ((dd->dim[0] == XX && j == YY) ||
- (dd->dim[0] == YY && j == XX))))
- {
- corner[j] += corner[i]*box[i][j]/box[i][i];
- }
+ corner[j] += corner[i]*box[i][j]/box[i][i];
}
}
if (c == 0)
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
+/* There's nothing special to do here if just masses are perturbed,
+ * but if either charge or type is perturbed then the implementation
+ * requires that B states are defined for both charge and type, and
+ * does not optimize for the cases where only one changes.
+ *
+ * The parameter vectors for B states are left undefined in atoms2md()
+ * when either FEP is inactive, or when there are no mass/charge/type
+ * perturbations. The parameter vectors for LJ-PME are likewise
+ * undefined when LJ-PME is not active. This works because
+ * bHaveChargeOrTypePerturbed handles the control flow. */
void ewald_LRcorrection(int start, int end,
t_commrec *cr, int thread, t_forcerec *fr,
real *chargeA, real *chargeB,
real *C6A, real *C6B,
real *sigmaA, real *sigmaB,
real *sigma3A, real *sigma3B,
+ gmx_bool bHaveChargeOrTypePerturbed,
gmx_bool calc_excl_corr,
t_blocka *excl, rvec x[],
matrix box, rvec mu_tot[],
tensor dxdf_q, dxdf_lj;
real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
real L1_q, L1_lj, dipole_coeff, qqA, qqB, qqL, vr0_q, vr0_lj = 0;
- gmx_bool bFreeEnergy = (chargeB != NULL);
gmx_bool bMolPBC = fr->bMolPBC;
gmx_bool bDoingLBRule = (fr->ljpme_combination_rule == eljpmeLB);
{
clear_mat(dxdf_lj);
}
- if ((calc_excl_corr || dipole_coeff != 0) && !bFreeEnergy)
+ if ((calc_excl_corr || dipole_coeff != 0) && !bHaveChargeOrTypePerturbed)
{
for (i = start; (i < end); i++)
{
/* Global corrections only on master process */
if (MASTER(cr) && thread == 0)
{
- for (q = 0; q < (bFreeEnergy ? 2 : 1); q++)
+ for (q = 0; q < (bHaveChargeOrTypePerturbed ? 2 : 1); q++)
{
if (calc_excl_corr)
{
}
}
}
- if (!bFreeEnergy)
+ if (!bHaveChargeOrTypePerturbed)
{
*Vcorr_q = Vdipole[0] - Vself_q[0] - Vexcl_q;
if (EVDW_PME(fr->vdwtype))
real *C6A, real *C6B,
real *sigmaA, real *sigmaB,
real *sigma3A, real *sigma3B,
+ gmx_bool bHaveChargeOrTypePerturbed,
gmx_bool calc_excl_corr,
t_blocka *excl, rvec x[],
matrix box, rvec mu_tot[],
"Take timings for this many steps in the benchmark runs" },
{ "-resetstep", FALSE, etINT, {&presteps},
"Let dlb equilibrate this many steps before timings are taken (reset cycle counters after this many steps)" },
- { "-simsteps", FALSE, etINT64, {&new_sim_nsteps},
+ { "-nsteps", FALSE, etINT64, {&new_sim_nsteps},
"If non-negative, perform this many steps in the real run (overwrites nsteps from [TT].tpr[tt], add [TT].cpt[tt] steps)" },
{ "-launch", FALSE, etBOOL, {&bLaunch},
"Launch the real simulation after optimization" },
CHECK(ir->bPeriodicMols);
if (ir->ePBC != epbcNONE)
{
- warning(wi, "Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
+ warning(wi, "Removing the rotation around the center of mass in a periodic system, this can lead to artifacts. Only use this on a single (cluster of) molecules. This cluster should not cross periodic boundaries.");
}
}
}
}
static void init_pull_coord(t_pull_coord *pcrd, int eGeom,
- const char *origin_buf, const char *vec_buf)
+ const char *origin_buf, const char *vec_buf,
+ warninp_t wi)
{
int m;
dvec origin, vec;
+ char buf[STRLEN];
string2dvec(origin_buf, origin);
if (pcrd->group[0] != 0 && dnorm(origin) > 0)
if (eGeom == epullgDIST)
{
+ if (pcrd->init < 0)
+ {
+ sprintf(buf, "The initial pull distance is negative with geometry %s, while a distance can not be negative. Use geometry %s instead.",
+ EPULLGEOM(eGeom), EPULLGEOM(epullgDIR));
+ warning_error(wi, buf);
+ }
+ /* TODO: With a positive init but a negative rate things could still
+ * go wrong, but it might be fine if you don't pull too far.
+ * We should give a warning or note when there is only one pull dim
+ * active, since that is usually the problematic case when you should
+ * be using direction. We will do this later, since an already planned
+ * generalization of the pull code makes pull dim available here.
+ */
+
clear_dvec(vec);
}
else
RTYPE(buf, pcrd->kB, pcrd->k);
/* Initialize the pull coordinate */
- init_pull_coord(pcrd, pull->eGeom, origin_buf, vec_buf);
+ init_pull_coord(pcrd, pull->eGeom, origin_buf, vec_buf, wi);
}
*ninp_p = ninp;
t_pull_group *pgrp0, *pgrp1;
t_pbc pbc;
int c, m;
- double t_start, tinvrate;
- real init;
+ double t_start;
+ real init = 0;
dvec dr;
- double dev;
+ double dev, value;
init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
md = init_mdatoms(NULL, mtop, ir->efep);
pull_calc_coms(NULL, pull, md, &pbc, t_start, x, NULL);
- fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
+ fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
for (c = 0; c < pull->ncoord; c++)
{
pcrd = &pull->coord[c];
fprintf(stderr, "%8d %8d %8d ",
pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1);
- init = pcrd->init;
- pcrd->init = 0;
-
- if (pcrd->rate == 0)
- {
- tinvrate = 0;
- }
- else
+ if (bStart)
{
- tinvrate = t_start/pcrd->rate;
+ init = pcrd->init;
+ pcrd->init = 0;
}
- get_pull_coord_distance(pull, c, &pbc, 0, dr, &dev);
- fprintf(stderr, " %6.3f ", dev);
+
+ get_pull_coord_distance(pull, c, &pbc, t_start, dr, &dev);
+ /* Calculate the value of the coordinate at t_start */
+ value = pcrd->init + t_start*pcrd->rate + dev;
+ fprintf(stderr, " %10.3f nm", value);
if (bStart)
{
- pcrd->init = init + dev - tinvrate;
- }
- else
- {
- pcrd->init = init;
+ pcrd->init = value + init;
}
- fprintf(stderr, " %6.3f\n", pcrd->init);
+ fprintf(stderr, " %10.3f nm\n", pcrd->init);
}
}
ewald_LRcorrection(fr->excl_load[t], fr->excl_load[t+1],
cr, t, fr,
- md->chargeA,
- md->nChargePerturbed ? md->chargeB : NULL,
- md->sqrt_c6A,
- md->nTypePerturbed ? md->sqrt_c6B : NULL,
- md->sigmaA,
- md->nTypePerturbed ? md->sigmaB : NULL,
- md->sigma3A,
- md->nTypePerturbed ? md->sigma3B : NULL,
+ md->chargeA, md->chargeB,
+ md->sqrt_c6A, md->sqrt_c6B,
+ md->sigmaA, md->sigmaB,
+ md->sigma3A, md->sigma3B,
+ md->nChargePerturbed || md->nTypePerturbed,
ir->cutoff_scheme != ecutsVERLET,
excl, x, bSB ? boxs : box, mu_tot,
ir->ewald_geometry,
srenew(md->massT, md->nalloc);
srenew(md->invmass, md->nalloc);
srenew(md->chargeA, md->nalloc);
+ srenew(md->typeA, md->nalloc);
+ if (md->nPerturbed)
+ {
+ srenew(md->chargeB, md->nalloc);
+ srenew(md->typeB, md->nalloc);
+ }
if (bLJPME)
{
srenew(md->sqrt_c6A, md->nalloc);
srenew(md->sigmaA, md->nalloc);
srenew(md->sigma3A, md->nalloc);
- }
- if (md->nPerturbed)
- {
- srenew(md->chargeB, md->nalloc);
- if (bLJPME)
+ if (md->nPerturbed)
{
srenew(md->sqrt_c6B, md->nalloc);
srenew(md->sigmaB, md->nalloc);
srenew(md->sigma3B, md->nalloc);
}
}
- srenew(md->typeA, md->nalloc);
- if (md->nPerturbed)
- {
- srenew(md->typeB, md->nalloc);
- }
srenew(md->ptype, md->nalloc);
if (opts->ngtc > 1)
{
char stmp[STRLEN];
ir->nsteps = nsteps_cmdline;
- if (EI_DYNAMICS(ir->eI))
+ if (EI_DYNAMICS(ir->eI) && nsteps_cmdline != -1)
{
- sprintf(stmp, "Overriding nsteps with value passed on the command line: %s steps, %.3f ps",
+ sprintf(stmp, "Overriding nsteps with value passed on the command line: %s steps, %.3g ps",
gmx_step_str(nsteps_cmdline, sbuf),
- nsteps_cmdline*ir->delta_t);
+ fabs(nsteps_cmdline*ir->delta_t));
}
else
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff