Free energy calculation as a function of temperature over phase transition

Question

I've run a molecular dynamics simulation (NVT) of a simple Lennard-Jones system of a few particles and I'm trying to calculate the Helmholtz free energy $\Delta F$ as a function of temperature over the phase transition. To properly sample the potential energy histograms I've run a replica-exchange MD simulation of 9 replicas, with approximately equal overlap for a series of temperatures of the bound system up to the unbound (melted) system.

However, as for the difference in free energy $\Delta F$ calculation I'm quite confused. I'm thinking of first trying with free energy perturbation (FEP) since this should be simplest and then attempting e.g. WHAM. From what I've read, FEP relies on some coupling parameter often denoted $\lambda$, two forms of FEP equations (Zwanzig equation) for the free energy are:

\begin{equation} \Delta F = F(\lambda = 1) - F(\lambda = 0) = -\frac{1}{\beta} \sum_{i=1}^{n} \ln{\langle e^{-\beta (U_{\lambda{i+1}} - U_{\lambda_{i}})} \rangle_{\lambda_{i}}} \end{equation}

called "windowing", where $n$ is the number of $\lambda$s sampled or intervals from 0 to 1, $U$ is the potential energy and $\beta$ the inverse temperature. The "integration" scheme: \begin{equation} \Delta F = F(\lambda = 1) - F(\lambda = 0) = \sum_{\lambda = 0}^{\lambda = 1} \left\langle \frac{dU_{\lambda}}{d\lambda} \right\rangle_{\lambda} \Delta \lambda = \sum_{\lambda = 0}^{\lambda = 1} \langle U_f - U_i \rangle_{\lambda} \Delta \lambda \end{equation}

$U_f, U_i$ are then the potential energy at the next $\lambda$ and the current, respectively? Similar to the summing above, in windowing.

I can sort of intuitively understand the process for when $\lambda$ represents how much one potential is "turned" on. Like when one gradually (by increasing $\lambda$) inserts a molecule in a solution and one then samples the energy during that process for different "amounts" of $\lambda$, this is the type of scenario I've read about in articles.

The problem I have is that I don't know exactly what my $\lambda$ is. Since I've only sampled over a range of temperatures and while this corresponds to my particles being closer or farther from each other my potential is the same over the whole simulation, only the temperature changes. In some sense then I believe my situation is simpler than the other and my naive approach was to simply take $\lambda$ to be the temperature $T$ in my case. And, for the first equation I took $\beta$ to be the temperature at $\lambda_i$ i.e. $T_i$ since it is here the average is supposed to be taken.

The two figures below are the results of my attempt at calculating $\Delta F$ but I'm quite sure neither is correct since they should be concave down ($\frac{d^2F}{dT^2} = -C_v/T < 0$ since $ C_v > 0$).

Using the first equation (windowing) my process was (in reduced LJ units):

• Calculate the average potential energy difference of each neighbouring temperature pair $T_{i+1}, T_i$ as $T_i \ln \langle exp(\frac{1}{T_i}(U_{T_i+1} - U_{T_i}))\rangle$ from my lowest temperature to my highest and sum these contributions. But I feel like I perhaps should take the difference in temperature as well here?

For the second equation, integration (less sure):

• Summed the pairwise average difference in potential energy, again between neighbouring temperatures and multiplied by the difference in temperature $\Delta T$ i.e. summing terms of $\langle U_{T_{i+1}} - U_{T_i} \rangle (T_{i+1} - T_i)$ from the lowest to the highest temperatures.

I might have misunderstood something essential when using replica-exchange but my impression is that it just improves sampling.

Edit: Since, (if correct), my perturbation is the temperature, meaning the kinetic energy changes with perturbation, shouldn't I consider the total energy $E$ not just potential $U$?

Edit 2: The figures are now updated, looking more like what I'd expect but I'm still not sure if they're correct/reasonable. Now scaled with the respective temperature at each point, included is also a comparison between full energy $E$ and potential only $U$. They're also both concave down now and look qualitatively similar. Though the integration seem to give half as large values. From what I can gather there should be a "break" at about the the phase change temperature as the phases here are in equilibrium and past this point one of the two curves takes over as the one lowest in free energy but I'm not quite sure I can see one.

Edit 3: From what I can read I've now realised, I believe, that integration (thermodynamic) is explicitly for systems in which the potential function changes with $\lambda$. Meaning it is not applicable for me in my case where only temperature changes? However, the first option (windom method) should still be viable.