How to Define the Heating Rate in the Heat Fix

I am using the heating fix

fix ID group-ID heat N eflux

My question is how to understand the heating rate eflux. I initially tried a value at order of magnitude 10^6 eV/ps, but the simulation very quickly crashed within a few time steps. Evidently I am not understanding how eflux is defined, so my heating rate is far too large. Smaller values have led to more stable performance.

This value of the heating rate I am using came from
Phillips et al, A two-temperature model of radiation damage in alpha-quartz, 2010,

where they used LAMMPS to model radiation with a heating rate of 205.76 kcal/(mole fs), which I converted to eV/(atom ps) and rescaled by 270 to account for the atoms in the heating control volume.

My questions are

(1) Does the user have to manually rescale the magnitude of eflux with the time step magnitude, or does LAMMPS automatically do this? I am using a time step of 1 fs = 0.001 ps.

(2) Does LAMMPS automatically rescale by N, the number of time steps between applying the heat fix? The User Guide implies this, but does not mention rescaling by the time step magnitude as well, i.e. Question 1.

(3) Does the user also have to manually rescale the magnitude of eflux with the number of atoms in the group where we apply the fix, or does LAMMPS also automatically do this? From my reading of the User Guide it seems that the user does have to rescale the value of eflux written in the input deck with the group size to make eflux extrinsic.

(4) As a more general question, how does LAMMPS use the value of eflux to determine how much energy goes into each atom (of a given group) whenever it needs to apply the heat fix?

(5) Are there any heuristic “rule-of-thumbs” to decide if an eflux magnitude is too large?

The relevant paragraph in the documentation has this text:

If eflux is a numeric constant or equal-style variable which evaluates to a scalar value, then eflux determines the change in aggregate energy of the entire group of atoms per unit time, e.g. in eV/ps for metal units. In this case it is an “extensive” quantity, meaning its magnitude should be scaled with the number of atoms in the group. Note that since eflux also has per-time units (i.e. it is a flux), this means that a larger value of N will add/subtract a larger amount of energy each time the fix is invoked.

No. note: eflux determines the change in aggregate energy of the entire group of atoms per unit time

Yes. note: Note that since eflux also has per-time units (i.e. it is a flux), this means that a larger value of N will add/subtract a larger amount of energy each time the fix is invoked.

You usually do not look at the heat added per atom, but the total heat added. So if you start from heat per atom, you need to multiply by the number of atoms.

For these details you need to look at the source code:

  • first LAMMPS determines the amount of heat (i.e. kinetic energy to be added)
  • then LAMMPS determines the factor by which velocities need to be scaled after subtracting the center of mass velocity sqrt((ke+heat-ke_com)/(ke-ke_com)
  • then LAMMPS applies the scale factor to the atom velocities w/o the center of mass velocity

I would look at the example(s) in the corresponding example folder and start with a similar magnitude.

I would start by looking very closely at the description in the paper you are referring to whether the heat added is the total heat or heat per atom.