## Bug #132

### Threshold for triggering "Large VCM" warnings arbitrary.

**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 over 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. over 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 13 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 13 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 about 13 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 about 13 years ago

Now fixed in head branch too.