[lammps-users] Langevin thermostat question

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps users

I’m trying to simulate a system interaction with an implicit solvent at different temperatures. I first equilibrate the system with npt at room temperature then scale the velocities to the desired temperature (50K more) and switch to nve with Langevin thermostat. However the damping constant in the Langevin thermostat turned out to be as high as 1500 (real units) after a calculation I made to properly describe my solvent.

I noticed that for example when I scale the velocities to 350K just before the Langevin run, the temperature drops to 250K at the start of the Langevin run and the thermostat does not increase it no matter how long I run (which is kind of unexpected to me since the thermostat should eventually reach the TStop at the end). I tested lower damping parameters and the thermostat acted as it should. I concluded that the damping parameter causes this drop in temperature.

My question is, to simulate the system at some desired temperature under this relatively high damping factor, should I scale the velocities to a higher value so that when the temperature drops, it reaches my desired value or there is another way ?

This technique works but I kind of find it a more of a workaround instead of being based on a scientific principle. Any suggestions to keep the temperature as desired when the solvent damping is high is highly appreciated

Thank you

2

Because of how fix langevin takes input, a higher “damping” number (really a damping time) corresponds to a lower friction coefficient. A damping time of 1.5ps corresponds to a friction coefficient of 0.67/ps. In some contexts this is indeed very small, but without knowing your specific system (protein? Ionic liquid? etc etc) it is impossible to give specific advice. But since Langevin dynamics controls temperature through friction, it makes sense that a too-low friction coefficient results in an inability to control temperature.

A rigorous way to make sure a Langevin simulation gives similar results to an explicit solvent simulation is to track how relevant degrees of freedom fluctuate in a benchmark explicit solvent simulation, and choose a damping coefficient such that those fluctuations “look similar” (e.g.have similar autocorrelation spectrum) in a Langevin simulation. I suspect that if you do this in your system you will find that a damping time much smaller than 1.5ps is appropriate.

Cheers,
Shern

3

It’s worth noting that the Stokes drag formula (used to relate diffusivity to viscosity in the documentation for fix viscous) is derived in the context of a particular hydrodynamics model. All the best finding a more suitable damping time for your Langevin simulations!

Cheers,
Shern

4

Thank you for the useful advices.

My system is a 2D graphene and I’m trying to investigate the presence of “air” as a solvent. That is the reason the friction is low. I did calculate the damping parameter by referring to the fix viscous command and relating the gamma parameter to the friction force in fix Langevin command.

5

Dear Lammps users

I’m trying to simulate a system interaction with an implicit solvent at different temperatures. I first equilibrate the system with npt at room temperature then scale the velocities to the desired temperature (50K more)

this is an odd procedure? why the scaling? why not equilibrated to the target temperature right away? after scaling the temperature, your system will no longer be equilibrated and thus will have to be re-equilibrated or your results will be tainted by the lack of equilibration. since you are running with fix npt that will also concern the volume, which will likely change with changing temperature.

and switch to nve with Langevin thermostat. However the damping constant in the Langevin thermostat turned out to be as high as 1500 (real units) after a calculation I made to properly describe my solvent.

1.5ps is a rather long time constant for the thermostat. how long is your production simulation time? how much time do you leave for re-requilibration?
after running with fix npt, how do you determine the (fixed) volume for your simulation box? do you just take the instantaneous value or an average over the cell lengths when the system has equilibrated?

here is another test you should make: equilibrate your system, but do NOT scale the temperature and then continue ONLY with fix nve. does that maintain the temperature well? what is the pressure? did the variable cell run properly relax it?

I noticed that for example when I scale the velocities to 350K just before the Langevin run, the temperature drops to 250K at the start of the Langevin run and the thermostat does not increase it no matter how long I run (which is kind of unexpected to me since the thermostat should eventually reach the TStop at the end). I tested lower damping parameters and the thermostat acted as it should. I concluded that the damping parameter causes this drop in temperature.

My question is, to simulate the system at some desired temperature under this relatively high damping factor, should I scale the velocities to a higher value so that when the temperature drops, it reaches my desired value or there is another way ?

i don’t think you should do any scaling at all but rather equilibrate to the desired target properties. what is the point of equilibration if you then just push your system out of equilibrium?

This technique works but I kind of find it a more of a workaround instead of being based on a scientific principle. Any suggestions to keep the temperature as desired when the solvent damping is high is highly appreciated

I think the main problem is in your rather unscientific simulation protocol.

6

“air” is an odd choice for “solvent” on the time and length scales of atomic level MD simulations.
it is an even stranger choice for an implicit solvent since it will incur a rather large error from the implicit solvent approximation.

as an exercise, please do the following: assume that air is actually only nitrogen molecules and assume they behave like an ideal gas.

based on those assumptions you can compute how many nitrogen molecules would be in your simulation volume. I expect that you will be surprised by how few there are.

7

Dear Dr. Axel

Thank you for the useful comments. They are extremely helpful.

I was able to later anticipate some of the issues you pointed out such as the volume of the box. I now use a thermostat (Langevin) and a barostat (berendsen) with nve.

Regarding the scaling, I found this procedure in two publications reporting results comparable to experiments. Therefore I kind of disagree with your statement saying that the procedure is unscientific. In addition, I referred to a similar question in the manual where the exact solution I mentioned was proposed to overcome the drop in temperature due to the damping effects.

Scaling the velocities by an extra 50K does not really yield a far point from equilibrium. The Langevin thermostat should actually yield equilibrium during the first few ps (given an appropriate damping parameter) which is something I tested and confirmed by setting the damping to 100fs.

Nevertheless, using npt to reach the desired temperature is a more stable approach but I believe there is no difference when it comes to changes in temperature as low as 50K. In fact, 50K is a very normal fluctuation range during the production run.

Regarding the test you suggested, I performed it and the temperature dropped sharply from 300 to almost 10K. Does that mean that the system was not well equilibrate in NPT ?

Thank you for pointing out the air suitability as an implicit solvent. I will have to look more into this.

8

Dear Dr. Axel

Thank you for the useful comments. They are extremely helpful.

I was able to later anticipate some of the issues you pointed out such as the volume of the box. I now use a thermostat (Langevin) and a barostat (berendsen) with nve.

I don’t think that that will avoid the issue of instantaneous versus average box size.

Regarding the scaling, I found this procedure in two publications reporting results comparable to experiments. Therefore I kind of disagree with your statement saying that the procedure is unscientific.

These days, if something about MD simulations is published, it doesn’t automatically give it credibility. One has to look at the context. It is possible to publish MD simulation based research in journals where the reviewers have rather little experience with (and/or pay little attention to) statistical mechanics principles. In this specific case, the scaling not only makes little sense from that statistical mechanical point of view, it is contrary to common sense, too. So there would have to be a very good justification for the procedure.
Also, you have to take under consideration that you are misunderstanding the description.

In addition, I referred to a similar question in the manual where the exact solution I mentioned was proposed to overcome the drop in temperature due to the damping effects.

Pplease point out the exact location. Perhaps the text is misleading. As I mentioned above, context matters a lot in this case. Typical cases are when you are using constraints or rigid bodies or when you start from a minimized (and thus 0K) starting configuration and want to anticipate the amount of kinetic energy converted to potential energy during equilibration.

Scaling the velocities by an extra 50K does not really yield a far point from equilibrium.

again, this depends on the context and it is unnecessary. there has to be a good justification for it. as I mentioned in my previous email, you are just destroying most of the effort invested in equilibration.

The Langevin thermostat should actually yield equilibrium during the first few ps (given an appropriate damping parameter) which is something I tested and confirmed by setting the damping to 100fs.

Nevertheless, using npt to reach the desired temperature is a more stable approach but I believe there is no difference when it comes to changes in temperature as low as 50K. In fact, 50K is a very normal fluctuation range during the production run.

that is a misconception. while the instantaneous temperature may fluctuate, its fluctuation stems from the transfer of energy between kinetic and potential energy. the total energy should stay constant. with scaling you add to the total energy and that is different.

Regarding the test you suggested, I performed it and the temperature dropped sharply from 300 to almost 10K. Does that mean that the system was not well equilibrate in NPT ?

such a sharp drop indicates a serious problem. there is something very wrong with your input(s) and/or procedure. there may be a small change but the first set of energies should be the same as the last set of the previous simulation when you are doing a proper restart or continuation.

Thank you for pointing out the air suitability as an implicit solvent. I will have to look more into this.

please do the computation of the number of gas molecules for your simulation volume. it is worth the effort.