I tried to reproduce the harmonic oscillator described in Tuckerman’s book by adding a chain as instructed, but I still got standard error results and bimodal distribution.

You are using a custom pair style that is not part of LAMMPS, so I cannot comment on whether it is correctly implemented. But you can use bond style harmonic for this just as well.

This command will add the same velocity to both of your atoms, which will then proceed to move - as a pair - through space. For a harmonic oscillator either one of the atoms has to be immobile or both atoms have of be initialized to equal magnitude but opposite sign velocities.

Please note that this will use an internal temperature compute and temperature computes in LAMMPS subtract 3 degrees of freedom by default as required for fully periodic bulk systems due to their translation invariance. For this special case, you need to add those back using compute_modify.

The following should be closer to what you are looking for. Here is only one mobile atom bound to a second, immobile atom. It also outputs the temperature as computed internally by fix nvt and used for thermostatting and contains the command to re-add the 3 removed DOFs for translation invariance (otherwise the computed temperature would be returned as 0 since there are no DOFs left when translation invariance is applied).

units lj
dimension 2
boundary p p p
atom_style bond
atom_modify map array
region box block -2 2 -1 1 -1 1
create_box 1 box bond/types 1 extra/bond/per/atom 1
create_atoms 1 single -0.5 0.0 0.0
create_atoms 1 single 0.5 0.0 0.0
mass 1 1.0
create_bonds single/bond 1 1 2
pair_style none
bond_style harmonic
bond_coeff 1 20 1
group mobile id 1
velocity all set 0.0 0.0 0.0
velocity mobile set 0.5 0.0 0.0
fix flat all enforce2d
fix xonly all setforce NULL 0.0 0.0
variable x1 equal x[1]
variable x2 equal x[2]
variable x equal x[2]-x[1]
variable zv equal vx[2]+vx[1]
variable vx1 equal vx[1]
variable vx2 equal vx[2]
variable v equal vx[2]-vx[1]
variable T equal temp
fix move mobile nvt temp 1.0 1.0 1.0 tchain 4
timestep 0.001
compute_modify move_temp extra/dof -3
thermo_style custom v_x1 v_x2 v_x v_zv v_vx1 v_vx2 v_v temp c_move_temp
thermo 1000
run 10000000

Thank you very much! I am new to LAMMPS. However, I just ran this code and found that the distribution is still not the correct canonical distribution. What could be the problem?

I am not an expert in this and have not studied the book in question, so I don’t even know what you are comparing to what. I was just pointing out fundamental issues with your input where it differs from what you claim it does.

I sincerely appreciate your assistance! I am truly grateful for your guidance on the matter of LAMMPS subtracting 3 degrees of freedom by default, which has also resolved my confusion regarding the absence of temperature for individual particles. I am attempting to verify a harmonic oscillator using Nosé-Hoover method is just like microcanonical ensemble, while coupled to a Nosé-Hoover chain thermostat will be canonical ensemble.