Potential bug/unexpected interaction between LINCS constraints and too skew box

User discussions
1

GROMACS version:2024.3
GROMACS modification: Yes/No

Hi all,
I’m planning to run a MARTINI simulation of a polymer in a skewed triclinic box using anisotropic pressure coupling. To prepare, I equilibrated the system for over 100 ns using semiisotropic pressure coupling, which completed without any issues.

However, when switching to anisotropic pressure coupling, I consistently encounter “LINCS WARNING” messages, which eventually lead to the simulation crashing. At first, I assumed this might be due to inadequate equilibration, but I noticed that the LINCS failure is always preceded by a warning: Triclinic box is too skewed. Can not fix PBC.

From the GROMACS documentation and manual, I understood that skewed boxes should be handled by remapping coordinates using box vectors. However, after reviewing the source code, it appears that box remapping to handle skewness is only triggered on steps when the neighbor list is updated.

To test this, I set nstlist = 1 so that the neighbor list is rebuilt at every step. With this change, the simulation runs without errors and the LINCS warnings disappear. This suggests that there may be an interaction between the box geometry and LINCS constraints on steps between neighbor list updates when the skewed geometry is not corrected. I’m not deeply familiar with the internals of the constraint solver, so maybe someone more familiar might have more of an idea of whats going on.

Since updating the neighbor list every step has a significant performance cost, I’m wondering if this this behavior is expected or is it potentially a bug?

Thanks for any help

2

This must mean that your box deforms very fast. We have a margin of 0.4% before warnings start to appear. What is your nstlist? What is your tau_p? Is your compressibility value reasonable? Do you have an anisotropic system?

3

I was using standard values taken from the martini database. My nstlist was 20 - I now realise this is probably too large. tau_p I tried a range of values to see if this alleviated the issue from 5 to 100 ps. My compressibility was 3e-4. My system is anisotropic.

4

An nstlist value of 20 is rather small. I would hope things should work with tau_p=100 ps.

In general anistropic coupling can show rapid changes in box components, even more so with an anisotropic system. Do you have fully anisotropic coupling, also with scaling of off-diagonal box elements? And if so, why?

5

Yes I am scaling the off diagonal components. The reason I am using this coupling is because I am trying to model an anisotropic protein crystal

6

What compressibility value did you set for the off-diagonal components?