if (bEquiv)
{
/* set index for matching atom */
- noe_index[j] = groupnr;
+ noe_index[i] = groupnr;
/* skip matching atom */
i = j;
}
This is supposed to cover all CH3 groups and the like */
anmi = *atoms->atomname[index[i]];
anmil = strlen(anmi);
- bMatch = i < isize-3 && anmi[anmil-1] == '1';
+ bMatch = i <= isize-3 && anmi[anmil-1] == '1';
if (bMatch)
{
for (j = 1; j < 3; j++)
"equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
"a set of equivalent atoms specified as residue number and name and",
"atom name; e.g.:[PAR]",
- "[TT]3 SER HB1 3 SER HB2[tt][PAR]",
+ "[TT]HB* 3 SER HB1 3 SER HB2[tt][PAR]",
"Residue and atom names must exactly match those in the structure",
"file, including case. Specifying non-sequential atoms is undefined."
}
/*do a first step*/
- natom = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
+ natom = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
+ teller = 0;
do
{
rmsnow = rms_diff(isize, d, d_r);
fprintf(fp, "%g %g\n", t, rmsnow);
+ teller++;
}
while (read_next_x(oenv, status, &t, x, box));
fprintf(stderr, "\n");
gmx_ffclose(fp);
- teller = nframes_read(status);
-
close_trj(status);
calc_rms(isize, teller, dtot, dtot2, mean, &meanmax, rms, &rmsmax, rmsc, &rmscmax);
warning_note(wi, warn_buf);
}
- if (ir->coulombtype == eelPMESWITCH)
+ if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
{
if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
{
- sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
- eel_names[ir->coulombtype],
- ir->ewald_rtol);
+ real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
+ sprintf(warn_buf, "The switching range for should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
+ percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
+ warning(wi, warn_buf);
+ }
+ }
+
+ if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
+ {
+ if (ir->rvdw_switch == 0)
+ {
+ sprintf(warn_buf, "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones potential. This suggests it was not set in the mdp, which can lead to large energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw switching range.");
warning(wi, warn_buf);
}
}
return ener_npt;
}
-static real vrescale_gamdev(int ia,
+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;
double rnd[2];
- if (ia < 6)
+ assert(ia > 1);
+
+ do
{
do
{
- x = 1.0;
- for (j = 1; j <= ia; j++)
+ do
{
gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
- x *= rnd[0];
+ v1 = rnd[0];
+ v2 = 2.0*rnd[1] - 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;
}
- while (x == 0);
- x = -log(x);
- }
- else
- {
- do
- {
- do
- {
- do
- {
- gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
- v1 = rnd[0];
- v2 = 2.0*rnd[1] - 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;
- }
- while (x <= 0.0);
- e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
+ while (x <= 0.0);
- gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
- }
- while (rnd[0] > e);
+ 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;
}
return x*r;
}
-static real vrescale_sumnoises(int nn,
+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 r, gauss;
-
- r = 2.0*vrescale_gamdev(nn/2, step, count, seed1, seed2);
+ const real ndeg_tol = 0.0001;
+ real r;
- if (nn % 2 == 1)
+ if (nn < 2 + ndeg_tol)
{
- gauss = gaussian_count(step, count, seed1, seed2);
+ 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;
+ r += gauss*gauss;
+ }
+ }
+ else
+ {
+ /* 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,
+static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
gmx_int64_t step, gmx_int64_t seed)
{
/*
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff