Question about how LINCS calculates stress

GROMACS version: 2022.2
GROMACS modification: No
Dear Community,

I’m using GROMACS to simulate polymers and study the stress relaxation using Green-Kubo relationship. The polymer is merely bead-spring model. I compared the results between two systems, namely the bond between beads using either a harmonic spring of force constant 200kT/b^2 or LINCS constraint. If I only include bond potential ( or LINCS) without any other potential, I get the results as shown in Figure 1.

The stress produced from LINCS algorithm didn’t really match the theory but after a vertical shift it kinda overlaps with the theory curve and that produced by harmonic spring. However when I tried to simulate the system with bending and intermolecular potentials, I got completely different results on longer time scale. (Figure 2)

This quite puzzles me. My understanding is that either harmonic spinrg or LINCS constraint should only have impact on the short time scale when each single bond relaxes. Even though they don’t overlap at the short time scale they should coincide at longer time. However my results showed the opposite. So I’m wondering how LINCS calculates stress and why this affects its relaxtion a lot.

(P.S. G(t)=\frac{V}{3k_BT} \Sigma_{xy,xz,yz}\langle \sigma_i(t) \sigma_i(0) \rangle, but I normalize it by \nu k_B T = \frac{N_{\text{chain}} k_B T}{V}, so the curve in the plot are actually only virial autocorrelation \frac{G(t)}{\nu k_B T} =\frac{1}{3N_{\text{chain}} }\Sigma_{xy,xz,yz}\langle \Xi_i(t) \Xi_i(0) \rangle where \Xi is the virial. The simulation was using Brownian solver and NVT ensemble so the kinetic contribution to the stress was ignored.)

Thank you in advance!

What do you mean with “so the kinetic contribution to the stress was ignored”?

I think I agree with your expectations. The virial with LINCS could be sensitive to the accuracy of the constraints. Did you experiment with the LINCS parameters?

Dear Hess,

What I meant was that normaly stress (or pressure) also has the contribution from kinetic energy. But in Brownian dynamics the kinetic energy is pointless so it’s excluded.

I tested more conditions including LINCS order 4 or 8, harmonic springs with a force constant of 200 or 1000, and the results are shown below:

It still seems that LINCS will give different relaxation from the harmonic spring.

It looks to me that it’s not unlikely that the bond results will converge to LINCS when increasing the force constant.

Dear Hess,

After a few more tests, a stiffer spring did capture the behaviors that LINCS gave, even though this is a puzzling phenomenon since the bond stiffness normally shouldn’t affect the long-time relaxation. Thank you very much!

I guess it could, indirectly, affect it because it affects the stiffness of angles.