Small rigid molecules with constraints

User discussions
,
1

GROMACS version: 2023.1
GROMACS modification: No

I modelled a small rigid-body nitrate ion (or NO3+) with 6 constraints.
However, I observed weird and large oscillations during a simulation with LINCS (and many errors). It seems (at least visually) fine or OK with SHAKE.
Could you please let me know whether this is reasonable?
Is there any other better way to model such small rigid molecules?

Here is some more details.
I add the constraints to all N-O distances as well as O-O distances.
The total number of constraints is 6. The rotational and translational degree of freedom is 6. The sum is 12, being equal to the total number of degree of freedom of NO3. So the molecule should be rigid.

However, the N-O distance changes from 0.105 to 0.14 with LINCS.
For Shake, it changes from 0.1249 to 0.1254. (Expected at 0.125.)

My topology file reads,

[ moleculetype ]
;Name
NO3 2

[ atoms ]
;  nr   type    resnr     residue  atom cgnr      charge       mass
   1    N     1           NO3        N    1      0.8603       14.0070
   2    O     1           NO3        O    1     -0.6201       14.0070
   3    O     1           NO3        O    1     -0.6201       14.0070
   4    O     1           NO3        O    1     -0.6201       14.0070

[ constraints ]
;ai  aj funct     b0
  1   2     1     0.125
  1   3     1     0.125
  1   4     1     0.125
  2   3     1     0.21650635094
  3   4     1     0.21650635094
  4   2     1     0.21650635094

Here is the time evolution of N-O distance with LINCS (black) and SHAKE (yellow).
(The ordinate is in Angstrom.)
image

Parameters used are

continuation             = no
Shake-SOR                = no
shake-tol                = 0.00001
lincs-order              = 12
lincs-iter               = 8
lincs-warnangle          = 30
2

This system is likely not compatible with LINCS.

But why do you want to have a fully rigid molecule? Not only does this cause technical issues, but I don’t see an advantage over using harmonic bonds or at least non-constrained angles.

3

Thank you for reply.
Is there any plausible reason or typical explanation of this incompatibility?

Flexible models are also fine to me, but I tried to reproduce a rigid model used in a previous study.

4

For highly connected molecules the eigenvalues of the constraint connection matrix will be larger than 1, which is not supported by LINCS. In general one wants to avoids such situations as other constraint solvers would be slow.

5

Thank you for the comments.