Inconsistent temperatures/kinetic energies between ASE atoms and thermo.out file

Hi all, I'm hoping you know something I don't here. I'm running reverse NEMD calculations on Silicon and noticed some strange differences between temperatures and kinetic energies when output by GPUMD and when read back into ASE using calorine's gpumd.io module.

I tested this by initializing a box with a Maxwell-Boltzmann distribution, recording the energy and temp, then running for a very short time (0.01 fs) to see how much they changed. It seems like GPUMD reads the correct velocities and masses into the simulation, then calculates a wildly different temperature and kinetic energy than ASE does.

Here is the code that reproduces the issue.

System setup:

from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from calorine.gpumd.io import read_xyz, read_thermodynamic_data
from calorine.calculators import GPUNEP

atoms = read('equilibrated.traj')
MaxwellBoltzmannDistribution(atoms, temperature_K=10)
initial_atoms = atoms.copy()
# ----------------------------------------------------------------------------- #
# Equilibrate to temp
run_dir = 'rNEMD'
temp_K = 10
timestep = 0.01
dump_interval = 1

md_params = [
('dump_position', dump_interval),
('dump_exyz', [dump_interval, 1]),
('dump_velocity', dump_interval),
('dump_thermo', dump_interval),
 ('time_step', timestep),
 # ('ensemble', ['npt_ber', temp_K, temp_K, 200, 0, 100, 1000]), # [start_temp, end_temp, tau_t (coupling constant), pressure, bulk modulus, tau_p (coupling constant)]
# ('ensemble', ['nvt_ber', temp_K, temp_K, 100]),
('ensemble', 'nve'),
 ('run', 1)]

calc = GPUNEP('/home/dawson-smith/Desktop/Research/OtherSoftware/Interatomic_Potentials/Dong_Si_nep.txt', 
            command='/home/dawson-smith/Desktop/Research/OtherSoftware/GPUMD/src/gpumd',
            gpu_identifier_index=0,
            directory=run_dir,
            atoms = atoms)

# run simulation for steps_per_cycle
atoms = calc.run_custom_md(md_params, return_last_atoms=True)

The next block checks differences between the thermo.out file, the dump.xyz file (ported back to an ASE Atoms), and the difference in velocities before and after the simulation.

import calorine.gpumd.io as calio

# check temperature from thermo output
thermo = calio.read_thermo(run_dir + '/thermo.out', len(atoms))
print('Thermo file temperature:', thermo['temperature'].iloc[-1])
print('Thermo file kinetic energy:', thermo['kinetic_energy'].iloc[-1])

# check temperature from ase structure and momenta
dump = calio.read_xyz(run_dir + '/dump.xyz')
print('ASE Temperature:', dump.get_temperature())
print('ASE Kinetic energy:', dump.get_kinetic_energy())

# check velocities before and after 0.01 fs:
vels = pd.read_csv(run_dir + '/velocity.out', sep = ' ', header = None)
assert np.all(np.abs(vels - initial_atoms.get_velocities()) < 5e-5), 'velocities are not equal'
assert np.all(np.abs(initial_atoms.get_velocities() - dump.get_velocities()) < 5e-5), 'velocities are not equal'

And here are the outputs:

Thermo file temperature: 807.20246227

Thermo file kinetic energy: 0.1043391073175

ASE Temperature: 7.74986429289144

ASE Kinetic energy: 0.008013976881653978

The velocities from the dump.xyz file, velocity.out file, and initial ASE Atoms object are the same, so I don't think it's a unit issue. The masses and momenta are equal as well, so there shouldn't be a difference. The issue can be replicated with multiple ensembles (see commented lines in first code block) and multiple NEP potentials.

Do you all have any ideas for why the two outputs are so far apart? Thanks!