fix nve/fix heat with tip4p potential

Hi Lammps users,

I am simulating thermal transport through an organic-inorganic junction wherein I have a SAM (Self-Assembled monolayer) of alkanedithiol molecules between two Au leads with a layer of H2O at one interface between the thiol endgroup and the Au surface.

I first relaxed the initial structures using 'fix nvt' at 300 K. This works fine - temperature remains stable. The problem arises when I use 'fix nve' and subsequently 'fix heat' to generate a heat flux through the system. I add and remove heat (equal amount) at the penultimate layers of my Au leads (direction of heat flow is through the junction). The total energy monotonically drifts and does not remain constant at all (or oscillate about any average value). Even if I give the value of heat as 0.0, the same thing happens. So I looked at a few 'fix nve' relaxations I had performed earlier (but for much shorter time periods and with no 'fix heat'), I noticed the same behavior even though I hadn't caught it then since the deviations were small.

I am using TIP4P potential for H2O with long range Coulombic solver (PPPM). I am using SHAKE for constraining the bonds/angles for the H2O molecules and my time step is 0.5 fs (I have tried 1 fs as well) for a total time of 10 ns for the NEMD run. I have read the previous posts on the forum (for example: where similar problems were raised - I have tried the recommended solutions like adjusting time step, checking physics of system etc. but nothing seems to work. If I don't have the H2O layer (so only a Au-SAM-Au system), this NEMD setup works fine.

Note: For the attached files, my initial structure is quite relaxed so I've started directly out with a quick nvt at 300K and then the NEMD runs.

Attached files:

1. Input script -
2. Output file - log_au_sam_water_nemd_300K_1.lammps (I've cleaned this up a bit)
3. Structure file -

Shubhaditya Majumdar
PhD Student
Department Mechanical Engineering
Carnegie Mellon University (8.4 KB)

log_au_sam_water_300K_1.lammps (2.07 KB) (255 KB)

Paul may be able to comment on the use
of fix heat and whether there are issues
with using it with TIP4P and SHAKE.


Since you’ve seen that a heat value of 0.0 still shows energy drift, the problem seems independent of fix heat. Still, I’d recommend verifying this by removing fix heat completely, and doing a simple NVE run to check energy drift. From what you’re saying, I presume you’ll still see energy drift, even with fix heat completely removed.

One possible cause of small energy drift is insufficient accuracy with the long range Coulombic solver. To test this, try running with a high accuracy setting. Somewhere between 1e-6 and 1e-12 should be good. And you’ll probably want to try using Ewald instead of PPPM for this extreme high accuracy case.

You could also crank up the accuracy of the SHAKE solve.

And you might check your LJ cutoff. If it is relatively short, that could introduce some energy drift.

Even so, you’ll likely still get some small amount of energy drift. Then you’ll have to decide how small is small enough.

Comparing versus your NVT run isn’t a fair comparison since the thermostat masks the energy drift. I think there’s a way to check for non-conservation in the NVT case too, but it isn’t as obvious as in the NVE case where you simply make sure that KE+PE is a constant as the simulation progresses.


Thanks for the inputs. Increasing the accuracies of SHAKE and the solvers look like a good starting point. Let me see how they go.