Variable temperature during molecular dynamics simulation

[AugustinLu]

During NVT and NPT molecular dynamics simulations, the temperature of the thermostat is set by AtomMove.IonTemperature and remains constant. In some cases (ex: melt and quench), it might be useful to allow the thermostat's temperature to change. For instance, a linear evolution of the temperature for heating or cooling.

This could be done by adding a keyword AtomMove.IonTemperatureEnd, which would define the final target temperature.

1. By default, AtomMove.IonTemperatureEnd = AtomMove.IonTemperature.
2. If AtomMove.IonTemperatureEnd != AtomMove.IonTemperature, then the thermostat temperature evolves linearly from AtomMove.IonTemperature to AtomMove.IonTemperatureEnd.

[davidbowler]

This seems like a very sensible idea. As you say, a simple scaling would be easy to add. Is it possible that more complex temperature programmes might be required at some point? If so, then we should plan for this as well.

If you would like, please feel free to make a branch from develop to start adding this functionality.

[AugustinLu]

Thank you for your kind reply. One other possible time evolution of the temperature is an exponential decay following Newton's law of cooling: $T(t) = T_{environment} + (T0 - T{environment}) e ^{-kt}$ In that case, the constant $k$ should be carefully defined, and the temperature will tend to (but never reach) the target final temperature.

For the time being, I guess it is more reasonable to only use the simple scaling method (i.e. constant cooling/heating rate) for molecular dynamics simulations.

I am creating a new branch to start adding this functionality. How can we link a branch to this issue?

[JianboLin]

How can we link a branch to this issue?

You can directly create a branch from development button in the right side of this page, or you can create a branch first, then search the created branch from development button.

[AugustinLu]

Thank you for your reply.

So far, I do not have the rights to create a branch from this page.

The branch containing the proposed changes is located at: https://github.com/AugustinLu/CONQUEST-release/commit/91955c66a3042eaa750e3d2ec427710bd4ee5da6

[tsuyoshi38]

I would like to discuss the strategy for the implementation of this cooling or heating process. Sorry for my late reply, but I have checked the pull request #253 recently and I think we need some changes about it.

As far as I understand, the present scheme defines the cooling rate from the initial and final temperatures with the number of MD steps. I believe it is much better to set the cooling or heating rate explicitly, because we need to change Conquest_input for each job if it is set from the initial temperature.

Note that we assume AtomMove.NumSteps as # of MD steps for each job (not the total # of MD steps),

MDn_steps = fdf_integer('AtomMove.NumSteps', 100 )

and AtomMove.IonTemperature as the fixed Target temperature.

temp_ion = fdf_double ('AtomMove.IonTemperature',0.0_double)

In the code, temp_ion needs to be updated during the heating or cooling process, but I don’t think it appropriate to use this keyword as an initial temperature.

Anyway.. I think what we need to do are the followings. A) introducing new keywords for heating or cooling process and B) add one file (md.temperature ?) to print out final temperature

[tsuyoshi38]

A) For the new keywords, I think we need

• MD.VariableTemperature: showing that the MD is for cooling or heating process F(default) or T
• MD.VariableTemperatureMethod : at present only “linear” (default)
• MD.VariableTemperatureRate : cooling or heating rate in T/fs (since the unit of timestep is femto second.)
• MD.InitialTemperature : initial temperature in the process (K)
• MD.FinalTemperature : final temperature in the process (K)

Please suggest the better names if some of them sound inappropriate or weird.

A0) Variable Temperature or not, method

• MD.VariableTemperature: showing that the MD is for cooling or heating process F(default) or T
• MD.VariableTemperatureMethod : at present only “linear” (default)

It may not be necessary, but I think it is better to prepare the logical keyword MD.VariableTemperature, just for declaring that we want to perform MD for cooling or heating. (We still need to say that the MD is NVT or NPT).

At present, we only prepared the procedure to change the target temperature at a same rate. But, we may want to implement other ways to control the temperature in the future. So, it is probably better to prepare the keyword MD.VariableTemperatureMethod now, though we only have keyword linear at present.

A1) Cooling or Heating rate:

• MD.VariableTemperatureRate : cooling or heating rate in T/fs (since the unit of timestep is femto second.)

For this, we need some discussion and agreements. First, I assume negative values should be used for cooling processes. (and positive values for heating processes.) But, this is just the definition, so the opposite way is also possible.
We can check this definition if we have MD.InitialTemperature and MD.FinalTemperature .

Second, I think the unit should be in K/fs (femto second). The value 1.0 means 10^15 K/s This is because the unit of time step (AtomMove.Timestep) is femto second.

Please let me know if you are unhappy with these suggestions.

A2) initial and final temperatures, and job control

• MD.InitialTemperature : initial temperature in the process (K)

• MD.FinalTemperature : final temperature in the process (K)

  MD.InitialTemperature may not be necessary if we use MD.IonTemperature for this. But I am afraid it may cause misunderstanding. MD.FinalTemperature will be used to stop CONQUEST when the target temperature reaches this temperature. (To judge this, the sign of MD.VariableTemperatureRate can be used. I mean, we need to check the present temperature is lower or higher the target temperature for cooling or heating process.)

The present temperature can be calculated from the initial temperature (in the present job) and the MD step (in the present job) with the values of MD.VariableTemperatureRate and AtomMove.Timestep. The question is how to get "initial temperature" in the present run. (See the next.)

[tsuyoshi38]

B) add one file (md.temperature ?) to print out final temperature

The problem is we do not have a way to get the initial (target) temperature when we restart the CONQUEST MD run. From the atomic velocities read from 'md.checkpoints', we can calculate the "instantaneous" temperature, but this is not the final target temperature (th%T_ext) in the last run. I am not completely sure, but I have noticed that md.checkpoints does not include the information of target temperature.

There might be two (or three) options to solve this.

1. B1) add the information th%T_ext in md.checkpoints
2. B2) prepare another file like md.temp only in the case of MD calculations for cooling or heating process.
3. B3) put the information in InfoGlobal.i00.

I think B2) is better than other two options. For B1), the format of md.checkpoint is changed and we cannot use the old files any more, although the change is only for MD simulations of heating or cooling process. For B3), it is a little complicated. This file includes the information of atomic velocities which also exists in md.checkpoint. But this file is basically for generating the matrix elements (like L or K matrix) and the file may not be generated with some specific option. So, it is not a good idea to use this file to include the information of the final temperature.

We can also get the information of the initial target temperature from the history of MD simulations (md.stats). But I assume that time step can be different in different jobs. Considering such cases, I think B2) is the best.

[AugustinLu]

Thank you for your valuable comments.

In revision f527f8bc7d9913c728e7129f568a278ec78e0974, an updated implementation is proposed, including the additional keywords.

The main change is that everything now revolves around the temperature change rate, thus AtomMove.NumSteps no longer has any impact on the cooling/heating rate.

1. New keywords and variables

Support was added for the following keywords:

MD.VariableTemperature        T
MD.VariableTemperatureMethod  linear
MD.VariableTemperatureRate    -0.5
MD.InitialTemperature         2000
MD.FinalTemperature           300 

For now, we are considering that the temperature evolves linearly. At a later stage, an exponential decay could be implemented, for long simulations where we would like to mimic the real cooling of a system.

A positive (resp. negative) value for MD.VariableTemperatureRate means heating (resp. cooling).

<html xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">

**Type** | **Internal variable** | **Location** -- | -- | -- Bool | flag_variable_temperature | Module **md_control** String | md_variable_temperature_method | Module **md_control** Float | md_variable_temperature_rate | Module **md_control** Float | md_initial_temperature | Module **md_control** Float | md_final_temperature | Module **md_control** Float | temp_change_step | _Internal_ **md_run**

temp_ion_end has been replaced by md_final_temperature.

2. Condition for stopping the simulation

Depending on the sign of md_variable_temperature_rate, the simulation will stop in function of the current thermostat target temperature $T_{ext}$ :

• Stop if $T_{ext}$ < md_final_temperature for md_variable_temperature_rate$<0$ (cooling)
• Stop if $T_{ext}$ > md_final_temperature for md_variable_temperature_rate$>0$ (heating)

3. Definition of target temperature

In the previous revision, the target temperature was defined at each step by an expression. Now it will be the previous target temperature + the change of temperature per step:

• $T{ext}[t]$ = $T{ext}[t-1] +$temp_change_step

4. Tested configurations

• [x] NVE (no effect)
• [x] NVT-NHC
• [x] NVT-SVR
• [x] NPT-NHC-PR
• [x] NPT-SVR-PR
• [ ] Others (Berendsen, MTTK)

5. Restart

I would go for the B2 approach, adding an additional file md.temperature. This ensures that the format of the existing files remains unchanged.

6. Planned additions and questions:

• Check that the initial and final temperatures are consistent with the temperature change rate (i.e. $T{final}-T{initial}$ must have the same sign as md_variable_temperature_rate)
• Do we need to add conditions for the thermostat changes to not be applied if the simulation is NVE ? It basically makes no change since it is not used.
• Is it better to set the current target temperature (at time step $t$) as a function of the previous time step $t-1$ or as a general expression?

[tsuyoshi38]

I have started testing your code (sorry for my slow response), and have had the problems related to these questions.

Check that the initial and final temperatures are consistent with the temperature change rate (i.e. must have the same sign as md_variable_temperature_rate)

I think the code should check this consistency. In addition, it should be warned if AtomMove.IonTemperature is set (not zero) with MD.VariableTemperature being set as true.

By the way, let me confirm one point. I assume that MD.InitialTemperature is not used at restarted jobs, even when it is set in Conquest_input. Am I correct? (the value read from md.temperature is used instead.)

Do we need to add conditions for the thermostat changes to not be applied if the simulation is NVE ? It basically makes no change since it is not used.

I think that CONQUEST should be stopped with appropriate messages in such cases. (I mean, when NVE is selected and MD.VariableTemperature is set as true.)

Is it better to set the current target temperature (at time step ) as a function of the previous time step or as a general expression?

It is probably better to use a general expression, but please consider that the cooling or heating rate may be different from previous jobs. (Does it answer your question ?)

[AugustinLu]

Thank you for your reply.

I think the code should check this consistency. In addition, it should be warned if AtomMove.IonTemperature is set (not zero) with MD.VariableTemperature being set as true.

The code will now abort if the temperature change rate is inconsistent with the initial and final temperatures. A warning message has been added in the latter case.

By the way, let me confirm one point. I assume that MD.InitialTemperature is not used at restarted jobs, even when it is set in Conquest_input. Am I correct? (the value read from md.temperature is used instead.)

Yes. In the case of a restart, even if MD.InitialTemperature is defined in the input file, its value will be discarded and replaced by value read from md.temperature. Also, I added a condition for CONQUEST to read md.temperature only if MD.VariableTemperature is true.

I think that CONQUEST should be stopped with appropriate messages in such cases. (I mean, when NVE is selected and MD.VariableTemperature is set as true.)

From the latest commit (5f6366b), CONQUEST will abort with an error message in these cases.

It is probably better to use a general expression, but please consider that the cooling or heating rate may be different from previous jobs. (Does it answer your question ?)

Indeed, with a cooling rate that may change at each restart, it should be better to change the current target temperature by the temperature step. In that sense, the current implementation should be satisfactory.

The changes discussed above are reflected by 5f6366b

[JianboLin]

Thank you for your commit.

I have two minority comments:

1. In the Check consistency, I think one more condition can be considered: When MD.VariableTemperature is True, but the difference of MD.InitialTemperature and MD.FinalTemperature is very small or same.

2. I think the code should check this consistency. In addition, it should be warned if AtomMove.IonTemperature is set (not zero) with MD.VariableTemperature being set as true.

temp_ion is set to be zero_double in fire_quenchMD and quenchMD, and 300.0_double for others, as default. How about we only warn with the message similar as: ignore the setting of AtomMove.IonTemperature when MD.VariableTemperature is True.

[AugustinLu]

Thank you for your comments.

1. In the Check consistency, I think one more condition can be considered: When MD.VariableTemperature is True, but the difference of MD.InitialTemperature and MD.FinalTemperature is very small or same.

What should be the behaviour, should this happen? A warning or an abort? From this hypothesis, the condition for completion would be met and the simulation would stop directly.

2. I think the code should check this consistency. In addition, it should be warned if AtomMove.IonTemperature is set (not zero) with MD.VariableTemperature being set as true.

temp_ion is set to be zero_double in fire_quenchMD and quenchMD, and 300.0_double for others, as default. How about we only warn with the message similar as: ignore the setting of AtomMove.IonTemperature when MD.VariableTemperature is True.

As this variable-temperature mode is for 'md' only, we may simply check if temp_ion is different from 300 and warn the user, should it be the case:

if (abs(temp_ion-300) > RD_ERR) then
    call cq_warn(sub_name,'AtomMove.IonTemperature is set and MD.VariableTemperature is true.')
end if

By the way, I wanted to return a longer message ("[...]is true, thus AtomMove.IonTemperature will be ignored."), but there seems to be a length limit for the warning message.

[JianboLin]

For 1, I think it only happens when the user forgot to set these two values. So, an abort may be considered.

For 2, ..., thus ignore AtomMove.IonTemperature may satisfy the length limit

[AugustinLu]

Please check if 35eaa5a56ee5d956300a8e134f018ee78565a18d solves these issues.

By setting the check for 1), the simulation will stop before initializing the rest, which would avoid wasting time.

[JianboLin]

The latest commit looks fine for me. Thank you for your modifications.

[davidbowler]

Closed by #253