I’m trying to simulate the interfacial tension between heptane and CO₂. For heptane, I’m using the TraPPE-UA force field where bonds are fixed. For CO₂, I’m using the flexible TraPPE model. So I have two types of molecules with different bond constraints in the same system.
Since I can’t use constraints=all-atom in the .mdp file (because that would also constrain CO₂, which should stay flexible), I defined the constraints for heptane directly in its topology file like this:
[ constraints ]
; ai aj funct length_A length_B
1 2 1 0.1540
2 3 1 0.1540
3 4 1 0.1540
4 5 1 0.1540
5 6 1 0.1540
6 7 1 0.1540
During the NPT equilibration, everything works fine. But when I extend the box in the z-direction (to introduce an interface) and run an NVT simulation to calculate the surface tension, I get this error:
Step 29202371, time 29202.4 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 2.006206, max 47.963226 (between atoms 3839 and 3840)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
3838 3837 90.9 0.1540 0.7131 0.1540
3839 3838 89.6 0.1540 3.3775 0.1540
3839 3840 90.0 0.1540 7.5403 0.1540
3837 3836 79.6 0.1540 0.0660 0.1540
Fatal error:
2 particles communicated to PME rank 30 are more than 2/3 times the cut-off
out of the domain decomposition cell of their charge group in dimension y.
This usually means that your system is not well equilibrated.
This mostly happens at higher pressures when the CO₂ mole fraction is 0.8 or higher, and sometimes even after long simulation times (like 29 ns). It’s been really frustrating — I’ve been stuck on this for days.
If anyone has advice or has faced something similar, I’d really appreciate any help or suggestions.