Paradoxical Density in OPLS-AA: Diesel mixture density is abnormally LOWER than its lightest compone

User discussions
, ,
1

GROMACS version:2022.03
GROMACS modification: Yes/No
Here post your question

Dear GROMACS Community,

I am currently working on simulating a complex diesel surrogate fuel system (a mixture of various hydrocarbons ranging from C9 to C24, including n-alkanes, iso-alkanes, naphthenes, and aromatics) and I have encountered a highly unphysical anomaly regarding the system density.

The Problem: When running NPT equilibrations using the OPLS-AA force field, the calculated macroscopic density of my entire diesel mixture turns out to be LOWER than the density of the lightest pure component in the system (pure C9 / nonane).

From a basic physical chemistry perspective, since the mixture contains a large proportion of heavy components and aromatics/naphthenes (which have much higher densities), the overall density should be significantly higher than pure C9.

What I have checked so far:

  1. Equilibration: The system is fully equilibrated. The density vs. time curve has completely plateaued.

  2. Box Visual Check: I visualized the trajectory in VMD. The simulation box is uniformly filled with molecules; there are no macroscopic vacuum bubbles or phase separations (droplets).

  3. Dispersion Correction: I have already enabled DispCorr = EnerPres in my NPT .mdp file.

Has anyone experienced a similar “density drop” issue with hydrocarbon mixtures in OPLS-AA? Could this be related to specific interaction parameters (like 1-4 scaling factors), topology generation for ring structures, or something else I might be missing?

Any suggestions or insights would be greatly appreciated! Thank you in advance!

Best regards, Scorpion

2

Have you tried running with different force fields.

3

Thank you for your reply!

To answer your question, I haven’t tried other force fields (like CHARMM or AMBER) for this specific system yet. Do you suspect that OPLS-AA might inherently struggle with this type of complex hydrocarbon mixture?

I want to share a very interesting finding that makes this even more baffling. I just ran validation tests on the pure substances (using the exact same .mdp settings and topology generation methods). For the pure components, the simulated densities actually came out slightly higher than their experimental values.

This makes the anomaly extremely puzzling to me: If the force field overestimates the density for the pure components (meaning A-A and B-B interactions are strong), why would the density drop so drastically when they are mixed?

Could this indicate a significant failure in the cross-interactions (A-B interactions, perhaps between the rigid aromatics and the flexible alkanes) generated by the mixing rules? Have you seen OPLS-AA fail like this in cross-interactions before?

I would really appreciate any insights you might have on this!

4

Where are you getting the parameters from ?

Have you looked at the following reference. It should not be too hard to try out charmm.

https://pubs.acs.org/doi/10.1021/ct200908r

5

Hi, thank you so much for the reply and for sharing this excellent reference!

To answer your first question: I obtained the initial OPLS-AA topologies (including charges and standard parameters) using the LigParGen server.

Regarding the paper you mentioned (Siu et al. on L-OPLS), I am actually very familiar with it! We have strictly adopted their approach to optimize the long-chain alkanes in our system. Specifically, we applied the L-OPLS modifications (halving the epsilon for methylene hydrogens) exclusively to the flexible aliphatic chains, while carefully retaining the standard OPLS parameters for the rigid naphthenic and aromatic rings to preserve their inherent densities.

Even with this L-OPLS correction, while the pure component densities are perfectly matching experimental values, the density of the mixture still drops significantly.

I noticed you suggested trying the CHARMM force field. I am very interested in this idea! Could you share a bit more about why you recommend CHARMM over OPLS-AA for this specific type of complex hydrocarbon mixture?

Do you suspect that CHARMM’s mixing rules (Lorentz-Berthelot) might handle the cross-interactions between rigid aromatics/naphthenes and flexible alkanes better than the geometric mean rule used in OPLS-AA?

I would really appreciate hearing your insights on this! Thanks again for your time!

6

Without your topology files/gro/mdp files it is difficult for me to judge. Are the residues in the gro files in the same order as the top file. Do you see any warnings from grommp ? do you see any lincs warnings due to any constraints , are there any constraints ? Have you inpsected the energies ? What do they look like ? very very negative ? There are so many ways to check your simulatons before coming to the conclusion that it is a force field issue or due 1-4 interactions ? if you suspect it is a force field issue a quick check with charmm/amber is not going to hurt. If you suspect 1-4 interactions, chamm-gui will also give you a top file so can verify if ligpargen and charmm provide the same output by inspecting the [pairs] section. GROMACS uses nrexcl 3 i.e, it excludes 1-4 interactions and re-includes it in the pairs section explicity. the coefficients of these are generated if you use gen-pairs yes option. Hope this helps

7

Hi, thank you so much for the incredibly detailed troubleshooting checklist! I really appreciate you taking the time to help me rule out those foundational issues.

To answer your questions: Yes, I strictly ensured that the residue orders match perfectly between the .gro and .top files. I also run a rigorous Energy Minimization (EM) before any NPT equilibration, so there were no LINCS warnings, no constraint violations, and the potential energies were stable and completely normal. For the 1-4 interactions, gen-pairs = yes and nrexcl = 3 were indeed correctly applied in the .mdp files.

To further isolate the issue, I even ran simulations of the pure components (e.g., pure n-decane, pure methyldecalin, and pure methyltetradecylbenzene) using OPLS-AA. The pure densities matched the experimental values beautifully.

But here is the most exciting update—I took your advice and tested the CHARMM force field (via CGenFF), and it was a complete game-changer!

I built a specific equimolar ternary mixture using those three representative molecules (flexible chain + rigid naphthenic ring + heavy alkylbenzene):

Under OPLS-AA, the macroscopic density of this mixture plummeted to ~741 kg/m³, exhibiting a massive and unnatural positive excess volume.

Under CHARMM36, the exact same ternary mixture yielded a stable density of 840 kg/m³!

Your intuition about testing CHARMM was absolutely spot on. It strongly indicates that the geometric mixing rule in OPLS-AA significantly underestimates the cross-interactions between large rigid rings and flexible aliphatic chains, preventing them from tightly packing. In contrast, CHARMM’s Lorentz-Berthelot mixing rules handle this “gap-filling” behavior flawlessly.

I am now officially migrating my entire system to the CHARMM force field. Thank you again for pointing me in the exact right direction—you literally saved my project!