Project

General

Profile

Bug #375

constraints = hangles is broken

Added by Peter Eastman almost 10 years ago. Updated almost 10 years ago.

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

Description

The "constraints=hangles" option is broken. The problem is that it generates too many constraints (more than the number of degrees of freedom available to be removed), leading to a singular constraint matrix. For example, consider the atoms

...-C-CH3

It generates a total of 10 distance constraints for those five atoms, one for each pair. Five atoms have 15 degrees of freedom, minus 10 constraints leaves only 5 residual degrees of freedom, which is too few. Nine constraints would be exactly right to make them move as a rigid body with six degrees of freedom. By adding one constraint too many, it makes the constraint matrix singular with the result that both LINCS and SHAKE fail to converge.

This can be fixed by reducing the number of constraints added by this option. There is no benefit to constraining all angles involving a hydrogen in any case. You really only care about ones of the form H-X-H or X-O-H. Those are the only ones needed to allow a 4-5 fs time step, and it avoids the problems described above while also not reducing flexibility as much.

The following code change implements this behavior. Starting at line 136 in topshake.c (the 4.0.5 version), it currently reads

if ((nshake == eshALLANGLES) ||
(count_hydrogens(info,3,ang->a) > 0)) {

Change that to

int numhydrogens = count_hydrogens(info,3,ang->a);
if ((nshake eshALLANGLES) ||
(numhydrogens > 1) ||
(numhydrogens 1 && toupper(**(info[ang->a[1]]))=='O')) {

History

#1 Updated by Erik Lindahl almost 10 years ago

Fixed+committed in master & release-4-0-patches.

Also available in: Atom PDF