LINCS going haywire for hydronium ion

User discussions
1

Houston, we have a problem I haven’t seen in over a decade, so my guess is that I am doing something seriously wrong. LINCS errors for the hydronium ion (perfectly pre-relaxed with flexible bonds) during EM. Here is my constraint-based topology:

[ moleculetype ]
; molname	nrexcl
H3O		2

[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
     1  opls_H3O   1    H3O     OW      1      -1.40
     2  opls_H3H   1    H3O    HW1      1       0.80
     3  opls_H3H   1    H3O    HW2      1       0.80
     4  opls_H3H   1    H3O    HW3      1       0.80

[ constraints ]
1  2  2       0.098
1  3  2       0.098
1  4  2       0.098
2  3  2       0.1619
3  4  2       0.1619
2  4  2       0.1619

[ exclusions ]
1	2	3	4
2	1	3	4
3	1	2	4
4	1	2	3

Any comments? Thank you!

2

I guess the eigenvalues of the constraint coupling matrix are larger than 1 which makes the LINCS matrix inversion through series expansion diverge.

SHAKE should work.

3

Thanks Berk. We found an intermediate solution (EM in flexible bond mode, MD with h-bonds constraints). Wouldn’t setting the constraint algorithm to SHAKE globally mess with water? I am guessing it shouldn’t, because they all use an explicit SETTLE directive, but just to make sure. Thanks again!