Tscale doesn't agree with reversible scaling method?

[ireaml]

Hi!

Thanks for developing such a useful package!

I have used the reversible scaling method to calculate free energies at high temperatures (switching from 100 K to 840 K). It works great and agrees with equilibrium thermodynamic integration. I wanted to switch to the tscale method because it would allow me to use the KOKKOS acceleration with lammps (not available for the pair hybrid/scaled). However, when comparing the results from ts and tscale, seems like there's something wrong with tscale (as shown below)?

I've tested using long equilibration and switching times for tscale (N_eq = 25000 and N_switch=250000 with dt=0.002 ps) but we get similar results as for the short switching time (shown in the figure above). My system is a Tellurium interstitial in CdTe but I've observed the same issue for bulk CdTe (with a different forcefield).

All settings between tscale and ts are the same. I've checked the equilibrated structure after the tscale switch from 100 K to 840 K and looks ok. The lammps input script I used for the tscale simulation is attached below (adapted from calphy.phase.temperature_scaling).

Wondering if you've seen this disagreement before or have an idea of what could be causing it?

Thank you!

Input lammps file for tscale simulation:

#Simulation variables -----------------------------#
variable         RANDOM equal 1
variable         iter equal v_RANDOM

# Simulation control parameters.
variable         n_eq equal 25000  # Equilibration time (calphy default: 25000 steps)
variable         n_switch equal 250000   # Switching time (calphy default: 50000 steps)
variable         t0 equal 100 
variable         tf equal 839  

variable         li equal 1
variable         lf equal v_t0/v_tf
variable         p equal 1.01325
variable         pf equal v_lf*v_p
variable         t_damp equal 0.1  # Nose-Hoover (calphy default: 0.1)
variable         p_damp equal 0.1  # Nose-Hoover (calphy default: 0.1)

# Initalizes the random number generator.
variable         rnd equal round(random(0,999,${RANDOM}))
#------------------------------------------------------------------------------#

#---------------------------- Atomic setup ------------------------------------#
units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map yes
newton on
timestep 0.002  # 2 fs

# Create atoms
read_data restart.data

# Define interatomic potential
pair_style mace no_domain_decomposition
pair_coeff * * MACE_model_swa.model-lammps.pt Cd Te
#------------------------------------------------------------------------------#


#----------------------------- Run simulation ---------------------------------#
# Remap box to get the correct pressure
# run 0
# change_box all x final 0.0 y final 0.0 z final 0.0 remap units box
# Performed in previous simulation (when calculating g at 100K)

# Set thermostat and run equilibrium.
# npt uses Nose/Hoover temperature thermostat and Nose/Hoover pressure barostat
fix               f1 all npt temp ${t0} ${t0} ${t_damp} iso ${p} ${p} ${p_damp}
print            "Equilibrating at ${t0} K"
# Create velocity and equilibrate
velocity          all create ${t0} ${rnd} mom yes rot yes dist gaussian
run               ${n_eq}
print             "Equilibration done"
unfix             f1

variable          step    equal step
variable          dU      equal pe/atoms
variable          lambda  equal ramp(${li},${lf})

# Forward integration
fix               f2 all npt temp ${t0} ${tf} ${t_damp} iso ${p} ${pf} ${p_damp}
fix               f3 all print 1 "${dU} $(press) $(vol) ${lambda}" screen no file ts.forward_${iter}.dat
print             "Starting forward integration"
run               ${n_switch}
print             "Finished forward integration"
unfix             f2
unfix             f3

# Equilibrate at final temperature
fix               f1 all npt temp ${tf} ${tf} ${t_damp} iso ${pf} ${pf} ${p_damp}
run               ${n_eq}
unfix             f1
# Check melting
dump              2 all custom 1 traj.temp.dat id type mass x y z vx vy vz
run               1
undump            2

# Reverse scaling
print             "Starting backward integration"
variable          lambda equal ramp(${lf},${li})
fix               f2 all npt temp ${tf} ${t0} ${t_damp} iso ${pf} ${p} ${p_damp}
fix               f3 all print 1 "${dU} $(press) $(vol) ${lambda}" screen no file ts.backward_${iter}.dat
run               ${n_switch}
print             "Finished backward integration"
unfix             f2
unfix             f3
print             "Random number used to initialise velocities ${rnd}"
print             "Finished simulation"
#------------------------------------------------------------------------------#