Effect of fix ave/chunk temp/bias options on thermo output from compute temp/com

[Ibrahim]

Dear All

I am computing the temperature of a nanoparticle in water using:

compute Tnp np temp/com

compute Twat water temp/com

thermo_style custom step c_Twat c_Tnp

and also computing radial temperature profiles using fix ave/chunk.

compute myChunksnp np chunk/atom bin/sphere 40 40 40 0.0 39 39 units box

I tested four approaches:

Case 1

fix profnp np ave/chunk 10 100 1000 myChunksnp temp bias Tnp file np_temp.profile

Case 2

fix profnp np ave/chunk 10 100 1000 myChunksnp temp file np_temp.profile

Case 3

compute TempChknp np temp/chunk myChunksnp temp com yes

fix profnp np ave/chunk 10 100 1000 myChunksnp temp bias TempChknp file np_temp.profile

Case 4

compute Tnp_profile np temp/com

fix profnp np ave/chunk 10 100 1000 myChunksnp temp bias Tnp_profile file np_temp.profile

When I checked c_Tnp from thermo_style, cases 1 and 4 give identical results and also identical profiles.

However, I observe differences in the reported c_Tnp values in cases 2 and 3.

The same behavior is observed for c_Twat. The corresponding commands for water are analogous to those for np and have been omitted for brevity.

Can fix ave/chunk, together with different temp/bias computes, affect the reported c_Tnp or c_Twat during a run, or should these values remain identical if the underlying trajectory is unchanged?

Also, should I rely on the c_Tnp and c_Twat values reported by thermo_style?

I am using LAMMPS version [29 Aug 2024].

Many thanks in advance,

Ibrahim

[akohlmey]

@Ibrahim It is very difficult to give any specific advice to your post and questions because there is not enough context and information about what you are simulating and how, as well as what specifically you want to learn from your simulation. The only observation that I have right now is that a 1 angstrom width for the binning shells is extremely small and thus it will be difficult to get meaningful statistics and results.

[Ibrahim]

Thanks for your reply.

I would like to determine both the nanoparticle temperature and the radial temperature distribution within the nanoparticle and the surrounding water.

For the radial temperature distribution, we are following the methodology reported in the article: http://dx.doi.org/10.1021/jp410054j.

After energy minimization and thermal equilibration using NVT followed by NPT, I apply a heat source to the nanoparticle and monitor the resulting heat transfer to the surrounding water.

Specifically, I observe that:

Cases 1 and 4 give identical c_Tnp temperatures and profiles.

Cases 2 and 3 produce different c_Tnp temperatures from Cases 1 and 4.

As an additional check, I repeated the calculations without the radial temperature profiling (chunk/atom, temp/chunk, and fix ave/chunk). In that case, the temperature reported by c_Tnp is the same as in Case 2.

My concern is whether the use of fix ave/chunk with different temp/bias options can modify or influence the temperature values reported by compute temp/com in thermo_style, or whether all computes should remain independent observers of the same trajectory.

If needed, I can provide a minimal input script that reproduces this behavior.

Many thanks in advance,

Ibrahim

[akohlmey]

@Ibrahim

If needed, I can provide a minimal input script that reproduces this behavior.

That would be helpful to better understand what you are doing.
You have to keep in mind that a person from the outside has no
insight into the details of your work, yet you are arguing as if this
was all very obvious.

With the reproducer, you should also check if there has been no
difference in the behavior of LAMMPS in the two years since your
LAMMPS version was released, i.e. make sure that you are not
observing an undesired behavior that was already changed.

[Ibrahim]

Thanks for your response.

I have attached the input script used for these calculations.

I will also try to install the latest LAMMPS version and check whether the same behavior persists.

#Input script:

units metal
atom_style full
boundary p p p
newton on
bond_style harmonic
angle_style harmonic
atom_modify map yes

read_data system_equi2.data

group np type 1
group water type 2 3

bond_coeff 1 19.5140 0.9572
angle_coeff 1 2.3850 104.52

write_dump all xyz structure_initial.xyz modify sort id element Ag H O

pair_style hybrid lj/cut/coul/long 10.0 eam/alloy
pair_coeff * * eam/alloy np.eam.alloy Ag NULL NULL
pair_coeff 3 3 lj/cut/coul/long 0.0044231 3.188
pair_coeff 2 2 lj/cut/coul/long 0.000 1.000
pair_coeff 2 3 lj/cut/coul/long 0.000 1.000
pair_coeff 1 3 lj/cut/coul/long 0.03906 2.916
pair_coeff 1 2 lj/cut/coul/long 0.000 1.000

kspace_style pppm 1.0e-4

##Minimization

#fix freeze_np np setforce 0.0 0.0 0.0

#thermo 100
#thermo_style custom step pe ke etotal fmax

#min_style fire
#min_modify dmax 0.01 vdfmax 0.0
#minimize 0 1.0e-2 10000000 100000
#write_dump all xyz structure_minim.xyz modify sort id element Ag H O
#write_data system_minim.data

##Solvent Relaxation
#reset_timestep 0
#velocity water create 300 4928459 mom yes rot yes dist gaussian

#thermo 100
#thermo_style custom step temp pe ke etotal

#fix relax_nvt water nvt temp 300 300 0.1

#run 500000
#write_dump all xyz structure_equi1.xyz modify sort id element Ag H O
#write_data system_equi1.data

#unfix relax_nvt

##Solvent Density Equilibration

#thermo 100
#thermo_style custom step temp press vol pe ke etotal

#fix relax_npt water npt temp 300 300 0.1 iso 1.0 1.0 1.0 dilate water
#run 500000

#write_dump all xyz structure_equi2.xyz modify sort id element Ag H O
#write_data system_equi2.data

#unfix relax_npt

#unfix freeze_np

##Heating

reset_timestep 0

velocity np create 300 4928459 mom yes rot yes dist gaussian

fix int_np np nve
fix int_wt water nve

##Tether (COM spring)

fix teth np spring tether 10.0 40 40 40 0.0

##OR recenter # another option for COM control

#fix lock_center np recenter 40.0 40.0 40.0 shift all units box

variable Power equal 900.0e-9
variable eflux equal v_Power*6.241509e6
fix optheat np heat 1 v_eflux

region r_in sphere 40 40 40 18.0 side out units box
region r_out sphere 40 40 40 39.0 units box
region r_shell intersect 2 r_out r_in
group outer dynamic water region r_shell every 100

compute Touter outer temp/com

compute_modify Touter dynamic/dof yes

fix bath outer langevin 300 300 0.1 4928459 zero yes
fix_modify bath temp Touter

compute Tnp np temp/com
compute Twat water temp/com

thermo_style custom step c_Twat c_Touter c_Tnp pe ke etotal
thermo 100

##RADIAL TEMPERATURE PROFILE

compute myChunksnp np chunk/atom bin/sphere 40 40 40 0.0 39 39 units box
compute myChunksWt water chunk/atom bin/sphere 40 40 40 0.0 39 39 units box

compute Tnp_profile np temp/com
compute Twat_profile water temp/com

fix profnp np ave/chunk 10 100 1000 myChunksnp temp bias Tnp_profile file solid_temp.profile
fix profWt water ave/chunk 10 100 1000 myChunksWt temp bias Twat_profile file liquid_temp.profile

##Run control

dump mydmp all custom 1000 dump_prod.lammpstrj id type x y z vx vy vz fx fy fz
dump_modify mydmp sort id

timestep 0.001
run 1200000

write_dump all xyz structure_heat.xyz modify sort id element Ag H O
write_data system_heat.data

Many thanks in adavnce,
Ibrahim

[Ibrahim]

I used the latest LAMMPS version (26.03) to check the effect of including the bias option in the ave/chunk command on the global c_Tnp value. Additionally, the radial profiles differ between the two versions. The plots of c_Tnp for cases 1, 2, and 5 with both versions are shown below:

[Ibrahim]

Many thanks in advance,
Ibrahim

[akohlmey]

Disclosure: I still don’t understand what is your specific problem and thus don’t feel able to give a meaningful answer on my own. I figured that it is perhaps my problem, and thus I have drafted this reply with the help of an AI assistant, which I had read through the relevant LAMMPS source code and documentation. I have reviewed the text and verified the quoted code paths myself. But as with all AI generated content, it can contain hallucinations, so I encourage you to verify each claim yourself.

@Ibrahim Thanks for posting the complete input deck and the systematic comparison; with that, the behavior can be explained exactly. The short answer: nothing is broken, and your computes are independent observers in the sense that none of them changes the definition of what another one reports. But the bias option of fix ave/chunk has an unavoidable side effect at the floating-point round-off level, and in MD that is enough to make two otherwise identical runs follow different - but equally valid - trajectories. Here is the mechanism, step by step.

1. compute temp/com in thermo never touches the trajectory

When thermo_style evaluates c_Tnp, the compute subtracts the group COM velocity from each atom’s velocity in a local temporary and accumulates the kinetic energy from that. It never writes to the global velocity array. So the thermo diagnostic itself is fully non-intrusive, and its temperature definition is identical in all of your cases.

2. fix ave/chunk ... temp bias <ID> temporarily modifies the velocities

The temp value with a bias is accumulated differently. In its end-of-step processing, fix ave/chunk calls the bias compute’s remove_bias_all() method, then sums m*v^2 over the now “thermal” velocities, then calls restore_bias_all() to put things back:

• for compute temp/com (your cases 1 and 4) that is v -= vbias followed later by v += vbias, applied in place to the stored velocity of every atom in the group;
• for compute temp/chunk ... com yes (your case 3) it is the same operation, except each atom gets its own chunk’s COM velocity subtracted and re-added.

In exact arithmetic this subtract/re-add round trip is the identity. In IEEE-754 double precision it is not: (v - b) + b is not guaranteed to reproduce v bit-for-bit. So on every sampling step (every Nevery = 10 steps in your setup) the velocities of the profiled atoms can be perturbed at the level of the last bit, i.e. ~1e-16 relative. Without the bias keyword (your case 2), fix ave/chunk only reads the velocities and the trajectory is untouched.

(This in-place remove/restore idiom is used in many places in LAMMPS; thermostats with a bias do the same, e.g. your fix bath ... langevin with fix_modify bath temp Touter. It is by design, so that all consumers can operate on the actual velocity array.)

3. Round-off plus chaos equals diverging trajectories

MD is chaotic: two trajectories that differ in the last bit diverge exponentially and are effectively decorrelated after a few hundred to a few thousand timesteps. From then on, the instantaneous c_Tnp values of the two runs differ by the normal thermal fluctuation amplitude, even though both runs are equally correct samples of the same ensemble. This is the same well-documented effect that makes LAMMPS runs differ bit-wise when changing the number of MPI ranks, the compiler, the FFT library, and so on (see the second bullet on the “Common issues” manual page linked below).

This sorts your observations exactly:

• Cases 1 and 4 identical: Tnp and Tnp_profile are two instances of the same compute style on the same group, so they produce bit-identical bias velocities. The round-off perturbations are identical, hence identical trajectories, thermo output, and profiles. Whether you reuse the thermo compute or define a dedicated one makes no difference.
• Case 2 identical to the run without any profiling: no bias means no write access to the velocities, so the trajectory is bit-identical to the unprofiled run. This is the cleanest confirmation of the mechanism, and you observed exactly that.
• Case 3 different from both: the per-chunk COM bias subtracts/restores different numbers than the global COM bias, so the round-off perturbation sequence differs, producing a third, distinct trajectory.
• Old vs. new LAMMPS version: any code change that reorders floating-point operations has the same consequence, so step-by-step trajectories from 29Aug2024 and the current version cannot be compared either.

4. What this means for your analysis

• Yes, you can rely on the c_Tnp / c_Twat thermo output - as statistical quantities. Compare time averages over the (quasi-)steady-state part of the run, with error bars, never instantaneous traces. All of your cases must agree within those error bars; if they do not, that would be a real problem worth reporting.
• The same applies to the version comparison of the radial profiles: with Nfreq = 1000 each profile snapshot is averaged over only 1 ps, which is far from converged - all the more so with 1 Angstrom shells (as noted before, the inner shells contain very few atoms). Average over a much longer steady-state window before concluding that the two versions disagree.
• The case 1/4 vs. case 3 profiles also differ for a physical reason, not just round-off: bias Tnp_profile removes one global COM velocity of the whole group from every atom, while temp/chunk ... com yes removes each shell’s own COM motion. These are two different definitions of “thermal temperature”; pick the one that matches the methodology of the paper you are following.
• One caveat with the per-chunk COM bias: fix ave/chunk does not subtract any degrees of freedom for a removed bias (see the DOF discussion and the adof/cdof keywords in its documentation). Removing 3 DOF per chunk while still normalizing by 3N per chunk noticeably and systematically lowers the reported temperature of thin shells containing only a few atoms.

Sources

• FixAveChunk::end_of_step(), TEMPERATURE branch with the remove_bias_all() / restore_bias_all() calls: src/fix_ave_chunk.cpp
• ComputeTempCOM::remove_bias_all() / restore_bias_all() (in-place global-COM subtract/re-add) vs. the local-temporary evaluation in compute_scalar(): src/compute_temp_com.cpp
• ComputeTempChunk::remove_bias_all() (in-place per-chunk-COM subtract/re-add; com yes marks the compute as usable as a bias): src/compute_temp_chunk.cpp
• Manual, “Common issues that are often regarded as bugs”, second bullet on round-off and trajectory divergence: 5.1. Common issues that are often regarded as bugs — LAMMPS documentation
• Manual, fix ave/chunk (the bias keyword, and the note that no DOF are subtracted for a bias; adof/cdof): fix ave/chunk command — LAMMPS documentation
• Manual, compute temp/com and compute temp/chunk: compute temp/com command — LAMMPS documentation , compute temp/chunk command — LAMMPS documentation

[Ibrahim]

Thanks for the detailed explanation, and sorry if my issue was not clear.

I had assumed that fix ave/chunk was primarily a post-processing command, so I expected it would not affect the trajectory or the temperature and energy reported by thermo_style.

I now understand that enabling the bias option affects the simulation at the floating-point round-off level and can introduce numerical differences that propagate into those quantities, and this behavior is expected and does not indicate a bug.

I will follow the recommendations regarding my analysis and focus on comparing statistically averaged quantities and properly converged profiles. Thanks.

Best regards,
Ibrahim

[srtee]

Refer to documentation:

The removal of the center-of-mass velocity by this fix is essentially computing the temperature after a “bias” has been removed from the velocity of the atoms. If this compute is used with a fix command that performs thermostatting then this bias will be subtracted from each atom, thermostatting of the remaining thermal velocity will be performed, and the bias will be added back in.

In MD temperature is calculated (roughly) as the sum of kinetic energy, but kinetic energy varies depending on your frame of reference (for example, we are orbiting the sun at a very high speed but we do not take that into account unless we need to go to space).

Compute temp/com subtracts the centre of mass velocity when calculating the kinetic energy of a system and thus it will give a slightly different temperature. Which answer is “correct” depends on whether you expect the group you are thermostatting to be moving for some reason. Some people simulate flowing fluids for which the mass flow kinetic energy should not be considered thermal, that gives them a good reason to use temp/com. It is impossible for me to know from afar what your situation is.