Project

General

Profile

Bug #132

Threshold for triggering "Large VCM" warnings arbitrary.

Added by John D. almost 13 years ago. Updated over 12 years ago.

Status:
Closed
Priority:
High
Assignee:
Erik Lindahl
Category:
mdrun
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

The code in mdlib/vcm.c appears to compare the kinetic energy against an arbitrary threshold of 1 (I'm
not certain what the units are) to generate a "Large VCM" warning:

for(m=0; (m<DIM); m++)
ekcm += sqr(vcm->group_v[g][m]);
ekcm *= 0.5*vcm->group_mass[g];
if ((ekcm > 1) || debug)
fprintf(fp,"Large VCM(group %s): %12.5f, %12.5f, %12.5f, ekin-cm: %12.5e\n",
vcm->group_name[g],vcm->group_v[g][XX],
vcm->group_v[g][YY],vcm->group_v[g][ZZ],ekcm);

Our group has experienced repeated large VCM warnings in various systems (including a well-behaved
box of Lennard-Jones atoms!). I suggest this cutoff either be computed in a system-dependent manner
(perhaps dependent on the size of the system?) or set to an extremely large number that is certain to be
an indication of catastrophic failure.

To compute a system-dependent threshold, you might first consider that the probability distribution of
the total kinetic energy for a system in the canonical ensemble has a simple analytic form. You can
choose a kinetic energy that is extremely unlikely to occur for the given number of atoms and atomic
masses based on this. If you don't have this worked out or readily available, I'd be happy to furnish the
equations.

Thanks!

- John Chodera

History

#1 Updated by Erik Lindahl almost 13 years ago

Hi John,

This is definitely a bug. No variable in Gromacs should ever be compared to an absolute number, with the
exception of the numerical under/overflow limits GMX_REAL_MIN/MAX.

IF you have the time and are interested in working out the equation it would be great. I agree it's
straightforward, but it will likely take a bit of testing to come up with reasonable default values for the
commonly used setups, and I have to other important bugs to fix this week :-)

#2 Updated by John D. almost 13 years ago

Hi Erik,

I'll see what I can come up with and get back to you in a few days. I think the equation for the kinetic
energy distribution is floating around in my thesis somewhere, and I'll try to come up with an easy way
to determine the "maximum likely kinetic energy threshold" from a given probability threshold.

Cheers,

- John

#3 Updated by Berk Hess over 12 years ago

I have also ran into this problem several times.
But the solution is not that easy as just comparing to a distrubition
(which should be a simply Gaussian distribution with sigma=1/2 kT,
as for any other degree of freedom).
The problem here is that it should be nearly zero, as we are removing vcm.
We only want to detect buildup of kinetic energy.

The build-up would depend on temperature, so comparing to kT
would be more reasonable.
But I think I have only experienced problems when using BD.
With plain MD or SD you should never get this warning,
unless your temperature is very high, or your system is crashing.

Berk.

#4 Updated by David van der Spoel over 12 years ago

Berk, do you mean comparing the EkinCM for the entire system to kT, i.e.
independent of the number of degrees of freedom? That should be easy enough to
implement.

#5 Updated by David van der Spoel over 12 years ago

I've fixed this in the 3.3 CVS by setting a max Tcm of 1K. I thought it wasn't worth the effort to add another parameter for this, and modify the tpr file. I haven't changed it in the head branch because this is substantially different code.

#6 Updated by David van der Spoel over 12 years ago

Now fixed in head branch too.

Also available in: Atom PDF