What is the most reliable way to calculate mixing enthalpy of binary solvent mixture using GROMACS?


I have simulated two-component systems of water and organic solvent with GROMACS. In a research paper from my field of study, authors calculate 'enthalpy of mixing' with the software in order to compare it with experimental mixing enthalplies. Under the Experimental Section of the paper, authors show that they do this according to this term:

enter image description here

A comment added in the text says the following:

"The internal energy, U, of the liquid mixture was obtained directly from the potential energy (see Tables S1 and S2) during the simulation, and yi is the molar fraction of component i."

I am interested in calculating same thermodynamical properties for my systems. I am aware that I can do that with g_energy function in GROMACS, or at least some parts of it (Pressure, Volume, Potential Energy).

My question is: Firstly, is this an appropriate way of calculating enthalpy of mixing in GROMACS via this equation? Secondly, how can one specify with the g_energy function which of the liquid's U(or potential energy), P and V is calculated if I initiate the function with .edr files, i.e. how can I choose to calculate i'th component's properties?

Reference for the image and quote:

Aguilera-Segura, S. M.; Di Renzo, F.; Mineva, T. Structures, Intermolecular Interactions, and Chemical Hardness of Binary Water-Organic Solvents: A Molecular Dynamics Study. J. Mol. Model. 2018, 24 (10), 292. https://doi.org/10.1007/s00894-018-3817-2.

  • 1
    $\begingroup$ BTW, the authors of that paper, according to your image, need an introductory course about scientific typography. $\endgroup$ – mhchem Jun 13 '19 at 14:46
  • $\begingroup$ mhchem, could you elaborate? Edit Do you mean there should be a space before gamma (or molar fraction in their case)? $\endgroup$ – Koryphae Jun 14 '19 at 5:16
  • $\begingroup$ @Koryphae There are several typographical flaws in that formula: 1. textual subscripts (liq, mix) should be upright; 2. minus sign doesn't look like one (vertically misaligned and has tight spacing on both sides); 3. use of asterisk * as multiplication sign; 4. inconsistent notations ("liq,1" but "liq2" without comma). Maybe @mhchem can find even more inaccuracies:) In defense of the authors, they use proper upright subscripts $(ΔH_\mathrm{mix})$ in the text, so it might as well be a journal editor to blame, but overall quality is indeed rather poor. $\endgroup$ – andselisk Jun 14 '19 at 5:32

GROMACS won't be able to calculate the enthalpy of mixing directly but using the energy function in GROMACS you will be able to get energy, pressure, and volume of the simulations.

Going through the paper that you posted it looks like the ran multiple simulation boxes for the different components and mixtures. They did this because you will not be able to get U,P, and V for individual components in GROMACS. So in order to calculate the enthalpy of mixing for a 2 component mixture they ran 3 simulations: one of the mixture, and one of each component by themselves.

Once you have the simulation boxes for the different components and mixtures it should be pretty straight forward to calculate the enthalpy of mixing.

Just a side note: I am always a little nervous calculating thermodynamic expressions using pressure from GROMACS because the pressure will fluctuate so wildly. So if you're going to make sure you're careful about your equilibration run and your NPT equilibration run as well.

Below are the GROMACS manual for energy which you would use to get U,P, and V from your simulation. The next is a tutorial for an NPT run from Virginia Tech and their tutorials are really well written and should help as a starting point and finally the last is a chat between a user and GROMACS developer discussing pressure fluctuations in NPT.




  • $\begingroup$ Thank you for your Answer, @dtg67. What is actually puzzling for me, is as far as I can understand, there is an option for g_energy to choose number of molecules (-nmol) for the calculations. If following the formula they suggest and you calculate mixing enthalpies for pure components and the mixture, shouldn't number of molecules in a system for pure components be exactly same i.e. if you have mixture of 2500 molecules then you calculate thermodynamic parameters for 2500 of both pure components as well and then make the calculations according to the formula of mixing enthalpy? $\endgroup$ – Koryphae Jun 14 '19 at 7:55
  • $\begingroup$ Authors have provided their system compositions in Supporting Information and they have all different sizes for all solvents they use. $\endgroup$ – Koryphae Jun 14 '19 at 7:56
  • $\begingroup$ Yes you should use -nmol but thermodynamically we don't need to have the same N in all the simulations. Looking at the SI I think that they ran all the simulations at the same volume which is why the Ns are different $\endgroup$ – dtg67 Jun 14 '19 at 15:06
  • $\begingroup$ Okay, yes I see it now, that the using the equation, it has molar fraction component to it. However, when I calculate potential energy in my pure ethanol system, without specifying with -nmol, I get value of -100 000 kJ/mol. The paper discussed has value of -40 kJ/mol. I guess this is when you specify for "-nmol 1"? But now question arises, what is appropriate -nmol for the mixture? Shouldn't they be comparable? $\endgroup$ – Koryphae Jun 15 '19 at 12:02
  • $\begingroup$ I can't comment on the differences between the result that your simulation has with literature value but all the -nmol flag does is divide energies by the number that you specify. $\endgroup$ – dtg67 Jun 15 '19 at 22:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.