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"
#------------------------------------------------------------------------------#

[srmnitc]

@ireaml Thanks for reporting. At first glance, everything looks ok with the script. I will run a few tests and get back to you. Could you please tell me how you did the integration after the run?

[ireaml]

Thanks for the quick response!

I just used an adapted version of the calphy.phase.integrate_reversible_scaling() setting scale_energy=False as exemplified in calphy.routines.routine_tscale](https://github.com/ICAMS/calphy/blob/main/calphy/routines.py#L376C5-L376C19) (and providing the corresponding values for the initial temperature, pressure and the free energy at the initial temperature. It's the same script I used for ts (but setting scale_energy=False for tscale).

import os
import numpy as np
from scipy.integrate import cumtrapz
import scipy.constants as const
from ase.io.lammpsdata import read_lammps_data

#Constants
h = const.physical_constants["Planck constant in eV/Hz"][0]
hbar = h/(2*np.pi)
kb = const.physical_constants["Boltzmann constant in eV/K"][0]
kbJ = const.physical_constants["Boltzmann constant"][0]
Na = const.physical_constants["Avogadro constant"][0]
eV2J = const.eV

def _integrate_rs(
    simfolder,
    f0,
    t,
    natoms,
    p=0,
    nsims=5,
    scale_energy=False,
    return_values=False,
):
    ws = []
    p = p/(10000*160.21766208)

    # The code below is adapted to read forward/backward files from different directory structures
    files = os.listdir(simfolder)
    fwd_files = [f for f in files if f.startswith("ts.forward_")]
    if fwd_files:
        print("Found forward files:", fwd_files)
        files_fwd = [f for f in os.listdir(simfolder) if f.startswith("ts.forward_")]
        files_bwd = [f for f in os.listdir(simfolder) if f.startswith("ts.backward_")]
        simfolders = [simfolder] * len(files_fwd)
    # else if we have Iter_i folders with each containing the .dat files
    elif os.path.exists(os.path.join(simfolder, "Iter_1", "ts.forward_1.dat")):
        folders = [f for f in os.listdir(simfolder) if f.startswith("Iter_")]
        simfolders = [os.path.join(simfolder, folder) for folder in folders]
        # Get forward file in each folder
        files_fwd, files_bwd = [], []
        for folder in simfolders:
            files_fwd.append([f for f in os.listdir(folder) if f.startswith("ts.forward_")][0])
            files_bwd.append([f for f in os.listdir(folder) if f.startswith("ts.backward_")][0])
    else:
        raise FileNotFoundError("No forward files found")
    # Get ist of numbers from filenames
    nsims = [int(f.split("forward_")[1].split(".")[0]) for f in files_fwd]
    print(f"Found {len(nsims)} iterations")
    for i in range(len(nsims)):
        simfolder = simfolders[i]
        fdx, fp, fvol, flambda = np.loadtxt(
            os.path.join(simfolder, files_fwd[i]), unpack=True, comments="#"
        )
        bdx, bp, bvol, blambda = np.loadtxt(
            os.path.join(simfolder, files_bwd[i]), unpack=True, comments="#"
        )

        if scale_energy:
            fdx /= flambda
            bdx /= blambda

        #add pressure contribution
        fvol = fvol/natoms
        bvol = bvol/natoms
        fdx = fdx + p*fvol
        bdx = bdx + p*bvol

        wf = cumtrapz(fdx, flambda, initial=0)
        wb = cumtrapz(bdx[::-1], blambda[::-1], initial=0)
        w = (wf + wb) / (2*flambda)
        ws.append(w)

    wmean = np.mean(ws, axis=0)
    werr = np.std(ws, axis=0)
    temp = t/flambda

    f = f0/flambda + 1.5*kb*temp*np.log(flambda) + wmean

    if not return_values:
        outfile = os.path.join(simfolder, "temperature_sweep.dat")
        np.savetxt(outfile, np.column_stack((temp, f, werr)))
    else:
        return (temp, f, were)

[srmnitc]

Thanks! I will take a look and get back to you. I should be able to take a look before the end of the week.