Fixing temperature using multiple thermostats in LAMMPS

Hi all,

I am working on implementing MD simulations with multiple thermostats to control the temperature of individual slabs. To achieve this, I divide the system into multiple slabs along the y-direction. This approach helps mitigate large local temperature deviations caused by latent heat release during free solidification of a pure metal. Below is the relevant portion of my script:

velocity  all create ${temp} ${rand}

label           loop
variable        i loop 28

variable        y1 equal -28+(${i}-1)*2
variable        y2 equal -28+(${i}+0)*2
variable        ss equal 1000+${i}
region          reg${i} block -3 3 ${y1} ${y2} -3 3 units box
group           grp${i} region reg${i}

fix             t${i} grp${i} langevin ${temp} ${temp} 0.1 ${ss}

next            i
jump            SELF loop

############################################################
fix             2 all npt temp ${temp} ${temp} 0.1 y 0 0 1.0
run             ${Nstep}

This script runs successfully. Is this script correct? Also, how to choose the best damping time for the Langevin thermostats?

Thanks!

No. Just look at your log file.

Can you please be more specific? The script runs successfully, so what details of the script do you think are incorrect? My main concern is that the damping time may not be appropriate, which may affect the atom dynamics.

Or could I use multiple NPT thermostats for these slabs? But LAMMPS does not seem to support this.

The fact that an input runs to completion is by no means an indication that it is correct. It can be performing a completely unphysical simulation and still run to completion. It only means that you have no syntax errors.

You are thermostatting atoms multiple times. That is completely unphysical and bogus.

There is lots of literature about thermostats, including langevin thermostats. You would need to study those to learn what is appropriate for your circumstances. There is no global “this is the right choice”.

That said, you entire setup of applyting multiple thermostats to multiple layers makes no sense with a thermostat like fix langevin, since that is applied to individual atoms anyway. So just one thermostat would suffice.

That makes even less sense.

For a good reason, you cannot have multiple fix commands changing the same box dimension at the same time. That is completely unphysical.

Thank you for the detailed explanation and prompt response.

My goal is to simulate a free solidification process at a temperature slightly below the melting point. I already have a script using a global NPT thermostat, where the system size is free to change in the solidification direction:

fix             2 all npt temp ${temp} ${temp} 0.1 y 0 0 1.0

This script works well and produces meaningful results. However, this method is not good enough because solidification releases latent heat, which may raise the local temperature near the solid-liquid interface well above the target value. There are several methods to address this issue, and the most effective one I found in the literature (see below) involves applying multiple thermostats to control the temperature of small regions within the system.

J. Monk, et al. Determination of the crystal-melt interface kinetic coefficient from molecular dynamics simulations." Modelling and Simulation in Materials Science and Engineering 18, no. 1 (2009): 015004.

Is there a convenient way in LAMMPS to implement this multiple-thermostat strategy?

I greatly appreciate your help.

Applying the additional langevin thermostats is in my personal opinion turning your simulation of a real physical behavior into a computer animation. The release of the heat is a physical effect, but using a thermostat to suppress that has no physical equivalent.

Feel free to ignore my opinion, but don’t expect any further comments from me.

Hi @MZhong,

Consider @akohlmey first answer and have a look at your log file. It will contain warnings that hopefully will make clear(er) to you why your setup is bogus.

Consider reading fix langevin and NPT doc pages as well as the How-to thermostat and barostat pages. They’ll show you how to use langevin correctly.

Hi @akohlmey,

I believe the issue has two layers.

The first concerns the physics. I agree that suppressing heat release is unphysical in principle. However, when the goal is to measure the kinetic coefficient—specifically, the slope of the solidification velocity versus temperature—it becomes important to accurately capture the local temperature near the solid–liquid interface. While using a global thermostat is more physically justified, the global temperature does not necessarily reflect the local conditions at the interface. Based on the methodology in the published paper I referenced, it seems reasonable to assume that we could suppress heat release without significantly affecting the atomic dynamics relevant to the kinetic coefficient.

The second layer is the technical challenge of implementing multiple thermostats in LAMMPS. I am not fixed on using the Langevin thermostat—that was just one approach I tested. I have tried several alternatives, but none worked. I initially thought this would be a straightforward task in LAMMPS, but it appears to be more complicated.

@Germain, regarding the log file: I couldn not find any warnings related to the Langevin thermostat in log.lammps; it appears to be working perfectly. My main concern was about selecting an appropriate damping time, which may affect the atom dynamics. I have also attached my LAMMPS script run.in for reference. I will certainly continue reviewing the documentation for more insights.

Thank you both for your time and support.

log.lammps (66.1 KB)
run.in (5.2 KB)

Indeed I was mistaking, the warning I was thinking about does not show up in this case.

However, the point remains, you are using two thermostats on your atoms which makes no physical sense. You should fix this before looking for the effect of damping parameters.

Thank you, @Germain, for the explanation. After doing some research online and asking ChatGPT, I’ve revised my script as follows:

velocity  all create ${temp} ${rand} 

############################################################
label           loop
variable        i loop 28

variable        y1 equal -28+(${i}-1)*2
variable        y2 equal -28+(${i}+0)*2
variable        ss equal 1000+${i}
region          reg${i} block -3 3 ${y1} ${y2} -3 3 units box
group           grp${i} region reg${i}

fix             t${i} grp${i} langevin ${temp} ${temp} 0.1 ${ss} zero yes

next            i
jump            SELF loop

############################################################
fix             2 all nph y 0 0 1.0

run             ${Nstep}

I have changed npt to nph to avoid specifying the temperature twice. To make sure that temperature control is applied to the entire system, I have also changed the box size to

region        box block -3 3 -28 28 -3 3

Does this change make sense?

Your input still has significant issues. It seems like you’ve skipped the LAMMPS 101 instructions and are trying to “vibe-code” your way through it while seeking verification here.

I doubt you’ll find anyone willing to double-check if ChatGPT’s advice makes sense under those circumstances. Long story short, ChatGPT is good at writing inputs that look reasonable but aren’t.

3 Likes

I see what’s happening here. It might seem like I simply fed my questions to ChatGPT and came here for verification. But no. What ChatGPT initially provided was not correct, even I can tell myself, so I spent hours testing, searching, and reading the official documentation.

As a LAMMPS beginner, I have found the documentation quite dense and difficult to read; the information I need is often buried in a sea of details. Despite its imperfections, ChatGPT has been helpful, because it provides useful information I want so that I can explore more in the documentation.

I understand that many people here are true experts, and I really appreciate the level of insight. That said, my expectation was simply to find clarification from someone with a bit more experience than I have. Based on a script shared by my colleague, which worked for a similar purpose, I initially thought I just needed to make a few small adjustments. I thought this aligned with the purpose of the “LAMMPS Beginners” tag. Perhaps this forum isn’t the best fit for questions in my case, but I do appreciate the responses.

1 Like

Hi @MZhong,

There are quite a few different issues to address here.

As background, we have had some recent unexpected experiences with ChatGPT, and in general the use of GenAI is still quite polarising in the academic community. I’m sorry that you experienced a bit of that backlash.

As far as I can see, in isolation, it is a reasonable experiment to use nph integration together with Langevin thermostats to avoid thermostatting the system twice. I can’t tell you whether the final result will be “correct” or “wrong” (more on that later).

I can still tell you that there’s a general problem with your “multiple Langevin” idea, in the details section below.

Re: Langevin thermostats

The Langevin thermostat is completely local – it acts on every single particle in the system purely based on the particle’s own velocity.

In your original citation, the authors used velocity rescaling thermostats which are “global” – they adjust particle velocities based on the total kinetic energy of the system. In any nonequilibrium situation a global thermostat would allow larger inhomogeneities in temperature, while limiting only the total KE. That’s why the authors, in order to have more aggressive and localised thermostatting, needed to define multiple regional velocity rescaling thermostats.

Such an approach should not be necessary for Langevin thermostatting. Particles should behave very similarly under a single “global” Langevin thermostat as compared to multiple “local” Langevin thermostats.

But my detailed talk highlights two issues with your request – not that it’s a bad question, more that (as you’ve suggested) the LAMMPS forum isn’t necessarily the best place to ask.

Firstly, the issue you are grappling with is primarily scientific. It is true that (1) LAMMPS has very many features (2) the manual is very long to read through and quite dense at points. However, the central issue – that you are not sure how to thermostat your system – is separate. The LAMMPS manual is meant to tell you what each thermostat does; it cannot tell you which thermostat is right for your system. You need to know what you want to accomplish to answer that question for yourself.

Secondly, the LAMMPS forum is primarily run by volunteers and we mostly have experience with developing and maintaining the LAMMPS source code, alongside whatever projects we use the code for. Thus, if the LAMMPS code does something and you expected it to do something else, this forum is a good place to ask. If the LAMMPS code does something and you aren’t sure whether that is scientifically appropriate or not, we will not be well-placed to answer your question. I happen to have some unofficial experience with nonequilibrium thermostatting – but I am not experienced at all at solid state physics, solid state MD, or phase transitions, plus simply not having as many years’ experience as other volunteers here. And even if I did have that experience, it would not be a substitute for the eventual process of peer review during publication.

So – I can appreciate your frustration. But I am not sure you will get clear answers here. You may be better off contacting the authors of the paper you cited. You should also read MD textbooks, especially Frenkel and Smit (and Evans and Morriss on nonequilibrium topics).

4 Likes

Hi @srtee,

Thank you for the detailed and patient explanation. I now realize there may have been some misunderstandings in our exchange. From a scientific perspective, I initially thought the issue was relatively straightforward—but if it’s more subtle than I assumed, that’s perfectly fine. After all, such complexity is part of the research process.

My main expectation when posting was to clarify whether my usage of LAMMPS commands was correct. Thanks to your explanation (which I genuinely appreciate), I now understand that the Langevin thermostat acts on every atom in the system and that applying multiple Langevin thermostats is unnecessary—and potentially incorrect. I also wasn’t aware initially that combining fix langevin with fix npt could be problematic, as it effectively sets the temperature twice, which may cause unintended behavior.

This is why, when I received feedbacks that there were issues with my script, I was confused—I already suspected there were problems, which is exactly why I sought help here.

It seems I did not express my goals clearly enough, which causes the misunderstanding. Regardless, I sincerely appreciate the responses and insights from everyone. I fully understand and respect that this community is made up of volunteers generously sharing their time and expertise, and I am very grateful for the support.

I think one issue is that in your scripts, the atoms in the group seem to have both the langevin and the npt/nph thermostats acting on them. You need to make sure that distinct fixes apply to distinct atoms.

So for eg., you would use Langevin (thermostat) + nve (time integration) on your group and nvt (thermostat & time integration) on the rest of your sample.

Edit: I suppose your second script does avoid thermostating, but is that the right ensemble you want to be simulating your bulk sample in?