I am currently performing a non-equilibrium molecular dynamics (NEMD) simulation using the fix ehex command in LAMMPS to introduce a temperature gradient in a water–ethanol mixture.
The system consists of 760 TIP4P water molecules and 40 ethanol molecules. The water molecules are constrained using fix shake to maintain rigid bond lengths.
The initial structure is provided in water-ethanol_vel.data, and the simulation is run using the input script soret.txt.
After the simulation, I extracted the z-direction temperature profiles for ethanol and water separately from the output dump file. The results are shown in the attached image.
As you can see, a clear temperature gradient is formed in the ethanol component, but the water temperature profile appears nearly flat, as if unaffected by fix ehex.
How can I ensure that the TIP4P water molecules also develop a temperature gradient under fix ehex?
Any advice or suggestions would be greatly appreciated.
I find it quite suspicious. As water and ethanol form a mixture, you should not see such difference between the two molecules’ temperature, shake or not. My guess is that the system size (4 nm) is too small for imposing such large temperature gradient.
Please also note that a timestep of 2 fs is likely too large for flexible ethanol molecules, and timesteps that are too larges are known to generate temperature issues.
Ethanol’s molecular mass is 2.5 times that of water – furthermore, it looks like your ethanol molecules are flexible while the water molecules are held rigid with SHAKE (giving it extra bond vibration DoFs).
Thus each ethanol molecule has a much higher heat capacity than a water molecule – that is, adding the same amount of thermalised random momentum to an ethanol molecule heats it much less than doing the same to a water molecule. By the ancient tradition of Argumentum Ad Handwavium, an NEMD simulation will take the path of least perturbation (compared to a reference equilibrium simulation) which in this case would be heating or cooling the ethanol molecules preferentially to the water molecules.
But I would also take @simongravelle’s observation very seriously – excessively large timesteps often result in equipartition problems (because the integrator can’t sample the fine vibrations needed to transfer heat between DoFs) and having ethanol molecules drastically un-thermalized relative to their surrounding water could certainly be considered poor equipartition.
In this situation I’d advise you to set up a reference system of the ethanol-water solution between two LJ walls, thermostatted to high and low temperatures respectively (for best “bang for buck” make it a “mirror system” – search the forums!). This should create a temperature gradient across the system, while each wall thermostat fix will have a scalar variable recording how much heat has been pumped in at the hot wall and pumped out at the cold wall. The advantage of this reference is that it is physically intuitive and you do not have to worry about the algorithmic details of something like fix ehex.
The system size, molecular model, and timestep in my simulation were chosen based on the 2005 paper by Carlos Nieto-Draghi et al., titled “Computing the Soret coefficient in aqueous mixtures using boundary driven nonequilibrium molecular dynamics.” In their study, both the temperature and concentration gradients were successfully formed. However, they employed the RNEMD (reverse NEMD or Müller-Plathe) algorithm rather than the ehex algorithm I’m using.
It’s possible that the simulation conditions in their study are not directly suitable for mine due to this difference in methodology.
As you suggested, I’m now planning to:
Increase the system size, particularly along the z-direction,
Reduce the timestep from 2.0 fs to improve accuracy and ensure better thermal equilibration.
Would 1.0 fs be sufficiently small for flexible ethanol molecules, or would you recommend going even lower?
Thanks again for pointing me in the right direction — it’s been very helpful in reassessing my simulation setup.
1 fs is definitely more reasonable, but it can sometimes still be too large: The physical way to choose the timestep, in the present situation, is to calculate the vibration frequency of the harmonic potential omega = sqrt(k / mu), where mu is the reduced mass mu = m1 m2/(m1+m2), and m1 m2 the masses of two atoms connected by a bond with constant k. With molecules like ethanol, the bonds involving a hydrogen atom are usually the ones vibrating at the highest frequency (thus the ones limiting the choice of timestep).
I’ve never done it, but keeping the C-H and O-H bond rigids could be a solution for increasing the timestep safely. I’ve seen it done with other code like GROMACS, I suppose its possible with LAMMPS and shake but I haven’t tried yet.
Thank you very much for your thorough and clear explanation. It was extremely helpful and deepened my understanding.
Your explanation that NEMD systems tend to follow the path of least perturbation — which in this case favors heat exchange with ethanol — was particularly convincing.
Previously, I attempted to apply a temperature gradient using fix langevin. While the temperature gradient was successfully formed, the Soret effect — i.e., the concentration gradient of ethanol — did not appear. The input file I used in that attempt was langevin.txt.
In the paper “Soret and mass diffusion measurements and molecular dynamics simulations of n-pentane–n-decane mixtures” (Andrea Perronace et al., 2001), the authors noted that thermostats and shake are usually incompatible. To address this, they developed a method in which thermostatting and shake were applied alternately within each timestep until convergence was reached. However, such an implementation is quite technically challenging for me, and I suspect this limitation might have prevented the formation of a concentration gradient in my own simulation.
For this reason, I’m very interested in your suggestion of using a mirror system with LJ walls. May I ask:
Can the mirror system approach work effectively even when water molecules are constrained by fix shake?
Alternatively, if I separate water and ethanol into different groups and apply fix ehex independently to each, would that help in generating a temperature gradient in water as well?
Thank you again for your valuable insights — I’m learning a great deal from this discussion.
Thanks a lot for the explanation — the connection between vibrational frequency and timestep makes perfect sense now.
You’re right that C–H and O–H bonds are likely the limiting factors in ethanol due to their high-frequency vibrations. I hadn’t considered constraining only those bonds with fix shake, but I’ll definitely look into whether that’s feasible in LAMMPS.
For now, I’ll try a 1.0 fs timestep and check if it improves stability and energy distribution.
