Dear all,
I have been using the “fix heat” to compute the thermal conductivity of silicate glasses with the “buck/coul/long” potential. (by setting up heat source and sink). The temperature and total energy are found to increase slowly with increased steps. Further observation include:
-
Before using “fix heat”, equilibrium under the NVE ensemble is achieved.
-
Using the ReaxFF potential can achieve good convergence of temperature.
-
Using higher precision of Ewald summation (orginal value 1E-6) can slow down the temperature increase, but can not solve the problem.
So I am not sure whether the failure of convergence is due to the flying ice cube effect since energy is added/removed per MD step. Or it is due to something else. Please give me some suggestions if possible. Any comment or suggestion would be highly appreciated.
Flying ice cubes are unlikely with a solid. I assume you
mean you are using fix heat on 2 reservoirs to induce
a temperature gradient. So do your temperature profiles
achieve steady-state? What do you mean that the
temperature is rising, since it is different everywhere?
You could also try other thermal conductivity methods
mentioned in Section 6.20 of the manual.
Steve
Hi Steve,
Thanks for your reply.
Yes, steady-state temperature profiles are obtained and the calculated results are in good agreement with experimental data.
However, the running average of the system temperature is found to increase very slightly, so is the total energy. (averaging the values obtained from “thermo_style temp etotal”).
Since the simulations are performed under the NVE ensemble, and the energy added to the heat source is equal to that removed from the heat sink, I think the total energy should have converged. Although the energy increase can be very small (~1 eV for an 1 ns simulation), I am still confused with this observation because using the ReaxFF does not have this problem at all.
I think there are several issues that may be responsible for this.
-
Calculation of the Coulomb force (Ewald summation) may induce small errors.
-
“Fix heat” occasionally induces unphysically severe interactions between atoms via velocity scaling. (Some atoms may move too fast for one step)
But I have no confidence in these assumptions. Could you give me some suggestions?
Thanks again for your reply 
Hi Ye,
Since the simulations are performed under the NVE ensemble, and the
energy added to the heat source is equal to that removed from the heat
sink, I think the total energy should have converged. Although the
energy increase can be very small (~1 eV for an 1 ns simulation), I am
still confused with this observation because using the ReaxFF does not
have this problem at all.
In a recent paper (An enhanced version of the heat exchange algorithm with excellent energy conservation properties | The Journal of Chemical Physics | AIP Publishing) we showed that
/fix heat/, i.e. the HEX algorithm, does not conserve the energy very
well and proposed an enhanced method which does. I implemented the new
method in LAMMPS as /fix ehex
/(LAMMPS Molecular Dynamics Simulator). If the energy drift you
observed in your simulations is on the order of ~1% per ns, there is a
good chance that the drift is caused by higher-order truncation terms in
the time integration, as explained in the paper.
If you would like to test this, you only need to replace the /fix heat/
commands in your input file by /fix ehex /and change the arguments
slightly. Have a look at the examples HEAT/in.lj_ehex and
HEAT/in.spce_ehex or the doc pages for more information.
Best regards,
Peter
PS: /fix ehex /is part of the RIGID package. You might have to recompile
with "make yes-rigid".
I agree with Peter - if you are seeing small drift over ns timescales,
then it could be the method itself. Try ehex which is one of
the alternate methods I mentioned as discussed in Section 6.20 of the
manual.
Steve
Hi Steve and Peter,
Thanks you for the suggestion. I should have noticed the fix ehex command earlier.