The issue includes both the MD simulation and model sections.
MD Simulation Section:
In the first stage, 4000 argon atoms are equilibrated at 90K in the NPT ensemble. In the second stage, all argon atoms in the equilibrated liquid argon system are heated under the NVE ensemble, using the heating command “fix heat ar heat 1 eflux” (where eflux=1e-8). The units used in this simulation are in ‘real’ units, so the unit of eflux should be [kcal/mol/fs]. The system temperature T_MD at different time is obtained from the simulation results.
Model Section:
The relationship between the heat flux and system temperature in model (SI unit) is known to be：
Q [J/s]=m_Ar * c_pl * d(T_model)/dt
where Q is derived from eflux. Using this, a temperature vs. time curve, ‘T_model over time’, can be generated for the system.
The conversion relationship between Q and eflux is as follows:
Q [J/s] = eflux[kcal/mol/fs] * 4000/N_A [mol] * 4184 [J/kcal] * 10^15 [s/fs]
where m_Ar [kg]=0.3995[kg/mol]*4000/N_A ，c_pl=1122[J/(kg·K)].The value of c_pl is taken at 90K.
The temperature vs. time curves from the MD simulation and model are compared. As shown in the figure below, the yellow solid line represents the model’s temperature curve, and the blue solid line represents the MD simulation’s temperature curve. The yellow dashed line is the trendline for the model results, while the blue dashed line is the trendline for the MD simulation results. The two solid lines do not overlap as expected, and there is a significant difference in the slopes of their trendlines. Theoretical analysis suggests that the two curves should approximately overlap, but the plotted results do not align with this expectation.
I am considering whether there might be an issue with my understanding of the fix heat command, or if the conversion process from eflux to Q involves problems related to the real unit system.
Could you please help analyze where the issue might be?
Thanks for the interesting query. It appears that the LAMMPS temperature is not changing at all, within statistical variation. Instead of worrying about units and so on, I suggest you try to figure how to make the LAMMPS temperature increase, and then worry about comparing to the model.
Thank you for taking the time to review my query and for your thoughtful suggestion.
I would like to clarify that during the heating process, the temperature in my LAMMPS simulations does increase. When the eflux value is relatively large, the temperature rise is more noticeable (from 91K to 140K). Conversely, when the eflux is small, the temperature increase is much slower (from 91K to 93K), which may appear negligible within statistical variations.
These graphs below show the comparison of MD (blue line) and mode (yellow line) for liquid heating process, each with different ‘eflux’ (for MD) and corresponding Q (for model). The yellow dashed line is the fitted line of the yellow solid line, and the blue dashed line is the fitted line of the blue solid line.
Here is a simple example that confirms that fix heat is doing what is written in the documentation. This is based on the “melt” example with four modifications:
I increased the cutoff to from 2.5\sigma to 5\sigma to significantly improve energy conservation
I switched the thermo output to turn off normalization (so it more closely matches other units settings)
I used more conservative neighbor list update settings
I added a variable econs to the thermo output, that outputs the current total energy minus the energy that would be added through fix heat and thus should remain (mostly) constant.
If you check the output, the value of v_econs remains close to the value of Etotal at the end of the 2000 steps of equilibration. This means that fix heat is adding exactly the amount of kinetic energy that is expected, i.e. 50\epsilon per time unit.
units lj
atom_style atomic
lattice fcc 0.8442
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create 3.0 87287 loop geom
pair_style lj/cut 5.0
pair_coeff 1 1 1.0 1.0
neighbor 0.3 bin
neigh_modify every 2 delay 4 check yes
fix 1 all nve
thermo 50
thermo_modify norm no
run 2000
reset_timestep 0
fix 10 all heat 1 50
variable econs equal etotal-dt*step*50
thermo_style custom step temp ke pe etotal v_econs
thermo_modify norm no
run 10000
Thanks Axel! Also, from Dr. Li’s new results for larger eflux (1e-5, 1e-6), it is clear that dT/dt is proportional eflux, as it should be. Dr. Li’s previous results for lower eflux (< 1e-7) showed dT/dt to be a constant, which clearly has nothing to do with fix heat (i.e. if the fix heat command was completely removed, the same small temperature rise would likely be observed). The discrepancy between LAMMPS and model is not our responsibility to figure out, but it is likely due to some missing constant factor, such as the number of atoms.
