# Multiscale modeling of dimerization thermodynamics of formic acid

Dominika O Wasik, Silvia Lasala, Olivier Herbinet, Konstantin Samukov, Sofía Calero, Thijs J.H. Vlugt

► To cite this version:

Dominika O Wasik, Silvia Lasala, Olivier Herbinet, Konstantin Samukov, Sofía Calero, et al. Multiscale modeling of dimerization thermodynamics of formic acid. Fluid Phase Equilibria, 2025, 594, pp.114356. 10.1016/j.fluid.2025.114356. hal-05064689

HAL Id: hal-05064689

https://hal.science/hal-05064689v1

Submitted on 12 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Multiscale modeling of dimerization thermodynamics of formic acid

Dominika O. Wasik<sup>a,b ID</sup>, Silvia Lasala<sup>c,*</sup>, Olivier Herbinet<sup>c</sup>, Konstantin Samukov<sup>c</sup>, Sofía Calero<sup>a,b</sup>, Thijs J.H. Vlugt<sup>d,*</sup>

## ARTICLE INFO

### Keywords:

- Molecular simulations

- Monte Carlo

- Umbrella sampling

- Equations of state

- Quantum Mechanics

- Reactive mixtures

## ABSTRACT

Heat pumps, which recycle waste heat, are a promising technology for reducing CO₂ emissions. Efficiently using low-grade waste heat remains challenging due to the limitations of standard heat exchangers and the need for more effective working fluids. This work introduces a multi-scale methodology that combines force field-based Monte Carlo simulations, quantum mechanics, and equations of state to explore the potential of formic acid as a new reactive fluid in thermodynamic cycles. Formic acid exhibits dimerization behavior, forming cyclic dimers in the gas phase, which can enhance the thermodynamic efficiency of heat recovery systems. The dimerization reaction of formic acid is crucial because it integrates chemical energy into thermodynamic processes, potentially improving the performance of heat pumps and other energy systems. The study implements umbrella sampling in Monte Carlo simulations to compute the thermodynamic properties of HCOOH dimerization, including equilibrium constants, enthalpy, and entropy. Results from two different methods to study dimer formation, namely the dimer counter method and the potential of mean force method, show strong agreement with the enthalpy of dimerization of -60.46 J mol⁻¹ and -62.91 J mol⁻¹, and entropy of -137.36 J mol⁻¹K⁻¹ and -146.98 J mol⁻¹K⁻¹, respectively. A very good agreement of the Monte Carlo results with Quantum Mechanics and experimental data validates the accuracy of the simulations. For phase equilibrium properties, the Peng-Robinson equation of state, coupled with advanced mixing rules, was applied and compared to Monte Carlo simulations in the Gibbs ensemble. This approach enabled the determination of the Global Phase Equilibrium of the system, vaporization enthalpy, phase composition, vapor and liquid densities of the coexisting phases, and entropy as a function of temperature. The agreement between the thermodynamic model and Monte Carlo simulations confirms the reliability of the methodology in capturing the phase behavior of the system. The findings demonstrate a promising approach for discovering and characterizing new reactive fluids, contributing to more efficient and sustainable energy technologies.

## 1. Introduction

The reduction of greenhouse gas emissions has become urgent as the European Climate Law requires the European Union economy and society to achieve climate neutrality by 2050 [1]. The global energy-related CO₂ emissions were estimated at 37.4 billion tonnes (Gt) in 2023 [2], which indicates a need for further advancements in sustainable energy technologies. The total energy-related emissions would have increased three times more than the actual increase of ca. 900 Mt between 2019 and 2023, if not for the application of five key clean energy technologies: heat pumps, solar photovoltaics, electric cars, wind power, and nuclear power [2]. Heat pump technology which operates

on a thermodynamic cycle that absorbs heat from a source, upgrades, and transfers it to a sink using compression work, has the potential for waste heat recovery towards the decarbonization of EU industry [3–6]. The most promising sectors for heat pump application are non-metallic minerals, non-ferrous metals, paper, and food, while chemical industry shows much lower potential due to the predominant need for high-temperature heat, where heat pumps are less efficient [4]. To reduce carbon emissions and utilization of non-renewable energy sources, heat pumps are being explored as an alternative to heating technologies based on fossil fuels. Heat pumps supplied only ca. 7% of the global building heating demand in 2021 [7]. According to the International

<footnote><sup>a</sup> Materials Simulation and Modelling, Department of Applied Physics, Eindhoven University of Technology, 5600MB Eindhoven, The Netherlands</footnote>

<footnote><sup>b</sup> Eindhoven Institute for Renewable Energy Systems, Eindhoven University of Technology, PO Box 513, Eindhoven 5600 MB, The Netherlands</footnote>

<footnote><sup>c</sup> Université de Lorraine, CNRS, LRGP, F-54000 Nancy, France</footnote>

<footnote><sup>d</sup> Engineering Thermodynamics, Process & Energy Department, Faculty of Mechanical, Maritime and Materials Engineering, Delft University of Technology, Leegwaterstraat 39, Delft 2628CB, The Netherlands</footnote>

<footnote>* Corresponding authors.
E-mail addresses: silvia.lasala@univ-lorraine.fr (S. Lasala), t.j.h.vlugt@tudelft.nl (T.J.H. Vlugt).</footnote>
Energy Agency (IEA), to achieve net zero emissions by 2050, the number of installed heat pumps must increase from the current 180 million to 600 million by 2030 [7]. Several countries, including Canada and the United States, encourage residents to adopt heat pump technologies for their space and water heating needs through, i.e. subsidies [7]. The analysis of Kosmadakis et al. [4] showed that the total potential of industrial heat pumps in the EU is 28.37 TWh/year which corresponds to 1.5% of the total heat consumption. The waste heat required to be recovered and upgraded by the heat pump to cover this consumption is approximately 21 TWh/year, representing 7% of the total waste heat potential in EU industries. Despite being a significant energy resource, waste heat is currently insufficiently utilized. Forman et al. [8] reported that around 72% of the world's primary energy supply is lost through exhaust/effluents (52%) and other losses (20%). The electricity generation sector accounts for the largest portion of this heat loss at 39%, with 88% of it being extremely low-grade heat (below 373 K).

Although harnessing waste heat is appealing, the low-temperature difference between the heat source and the working fluid hinders effective heat transfer in standard heat exchangers, affecting the efficiency [9]. While many heat transfer enhancement solutions exist, these are generally designed for medium-to-high temperature ranges (> 500 K) [10], with limited options for low-temperature applications [11]. Heat recovery techniques are broadly classified into two categories: active methods, which involve mechanical or electrical systems requiring external energy inputs, and passive methods, which utilize natural processes or simple design features without external energy inputs. Passive techniques are considered more competitive than active methods for recovering very low-grade heat (below 373 K) cost-effectively. Therefore, there is significant motivation to re-evaluate a potentially cost-effective yet under-explored passive method: heat transfer intensification through gas dissociation [9]. Currently, reactive working fluids that undergo endothermic dissociation reactions when heated or compressed and exothermic recombination when cooled, instantly progressing towards chemical equilibrium, are investigated as an alternative to inert working fluids in thermodynamic cycles, due to their higher energy-efficiency potential [3,12–14]. As the temperature and pressure change during various unit operations, e.g. turbines, heat exchangers, pumps, and compressors, the fluid can undergo spontaneous, rapid, and reversible chemical reactions. These reactions incorporate chemical energy into the thermodynamic processes, significantly influencing the efficiency of the thermodynamic cycle and the design of components such as turbines [15]. Several studies have shown that the thermodynamic efficiency of supercritical (Brayton) and subcritical (Rankine) cycles using reactive working fluids, particularly N₂O₄, is higher than that of cycles using specific inert fluids such as water, ammonia, helium, or CO₂ [15–18]. The main disadvantages of Brayton cycles with ideal gas working fluids are the substantial compression work relative to the net output, and the highly sensitive performance to the efficiency of the turbomachinery and pressure losses [19]. In studies that consider fictive reactive fluids to examine how reaction stoichiometry, enthalpy, and entropy affect the performance of Brayton power cycles and heat pumps compared to pure fluids, it was shown that reactive working fluids can improve power plant efficiency by up to 30% [12] and more than double the coefficient of performance of heat pumps based on a gaseous inversed Brayton cycle [3]. The general reaction for fictive reactive fluids is:

$$ \left(\frac{n}{m}\right) A_m \rightleftharpoons A_n \qquad (1) $$

where n, and m specify the number of atoms of the two molecules.

Recent theoretical studies are focusing on the use and characterization of reactive fluids, including N₂O₄/NO₂, Al₂Cl₆/AlCl₃, and Al₂Br₆/AlBr₃ [15], to enhance the performance of unit operations and thermodynamic cycles. Among these, N₂O₄ is the most thoroughly studied, with its thermodynamic properties being determined through equations of state [13,20–24], direct numerical simulation [9], Monte Carlo simulations [13,25], and Quantum Mechanics [13]. In the study

of Zhang et al. [9] it was shown that the N₂O₄/NO₂ chemical reaction can enhance the heat transfer coefficient by up to 600% compared to non-reactive fluids. Despite this potential, the broader adoption of reactive working fluids in thermodynamic cycles is limited by the lack of similarly reversible fast chemical reactions, aside from N₂O₄. Addressing this challenge, the European Research Council funded project REACHER [26,27] that aims to discover new working fluids, characterize their thermodynamic and kinetic properties, optimize thermodynamic cycle architectures, and validate performance calculations through experimental micro-power plant testing.

The formic acid (HCOOH) reacting system has the potential as a new reactive fluid in thermodynamic cycles due to its strong dimerization behavior [28]. Formic acid is the simplest carboxylic acid, serving as a prototypical molecule and one of the most stable neutral complexes, with a dimerization energy of ca. 63 kJ mol⁻¹ [29]. Formic acid is a significant product in the chemical industry, with a production capacity reaching 800,000 tons annually as of 2017 [30,31]. The global market value of HCOOH is projected to grow from 1.8 billion dollars in 2023 to 2.8 billion dollars by 2033 [32]. Producing HCOOH offers an economically viable method for reducing CO₂ emissions through direct capture at the production source and subsequent conversion [33]. In 2021, HCOOH was utilized in various industries, including animal feed and grass silage (27%), leather tanning (22%), pharmaceuticals and food chemicals (14%), textiles (9%), natural rubber (7%), drilling fluids (4%), and other applications (17%) [34,35]. The predominant method for HCOOH production involves methanol carbonylation, resulting in a formate ester that is subsequently hydrolyzed to HCOOH [34,36–38]. Alternative methods for HCOOH production include the electrochemical reduction of carbon dioxide in aqueous electrolyte solutions [39–45] and the hydrogenation of CO₂ [31,46–48].

The formation of cyclic HCOOH dimers that occurs in the gas phase is significant to the human metabolism [49] and applications, e.g. the oxidation of unsaturated hydrocarbons [50] and atmospheric chemistry [51]. The graphical representation of the formic acid dimerization equilibrium is shown in Fig. 1.

The dimerization of formic acid (HCOOH) has been extensively studied both theoretically and experimentally [28,29,52–57]. While the studies by Jedlovszky and Turi [58] and Schnabel et al. [59] provide valuable insights into the modeling of formic acid systems, the focus of these works does not include the thermodynamic properties of HCOOH dimerization or the phase composition in terms of monomers and dimers. In this study, we focus on the thermodynamic properties essential for the design of heat pumps using reactive fluids like HCOOH, including vapor pressure, vaporization enthalpy, phase composition, and the densities of coexisting vapor and liquid phases. These properties, influenced by dimerization behavior, directly affect the heat transfer and energy storage capabilities of HCOOH, enhancing its potential as a working fluid in various industrial and engineering processes [28]. While HCOOH remains stable at room temperature, it exhibits thermal instability and decomposes when subjected to heat [60]. Despite extensive studies, there are still gaps in experimental data, particularly at varying temperature and pressure conditions, which makes the prediction of entropy and phase behavior challenging. Given the complex nature of HCOOH and the limited availability of experimental data, molecular simulations provide an alternative for determining these thermodynamic properties. Slavchov et al. [28] proposed a model to describe the temperature dependence of vapor pressure and heat of vaporization of associated liquids containing dimers and linear associates in the gas phase. The analytic generalization of the Clausius–Capeyron equation resulted in an enthalpy of dimerization of –58.5 kJ mol⁻¹ for formic acid. In a Monte Carlo simulation study by Turner et al. [52], the enthalpy of dimerization for formic acid in a carbon dioxide solvent was predicted to be –50.3 kJ mol⁻¹ in the gas phase. In a study by Chao and Zwolinski [60], several experimental values of the HCOOH dimerization enthalpy were reviewed. The average enthalpy of dimerization obtained from P-V-T methods is ca. –59 kJ mol⁻¹, while the
<figure><img src="image_4.png" /><figcaption>Fig. 1. Graphical representation of the formic acid dimerization equilibrium, given by a reaction: 2HCOOH $\rightleftharpoons$ (HCOOH)₂.</figcaption></figure>

average obtained from spectroscopic methods is ca. –57 kJ mol⁻¹. The chemicophysical properties of formic acid are influenced by its ability to form various types of hydrogen bonds (H-bonds) [53]. A formic acid molecule can act as a hydrogen acceptor via its two oxygen atoms and as a hydrogen donor through its hydroxy and formyl hydrogen atoms [53]. Consequently, four types of H-bonds are theoretically possible: O–H···Oc, O–H···Oh, C–H···Oc, and C–H···Oh, where Oc and Oh represent the carbonyl- and hydroxy oxygen atoms, respectively. This variety leads to different structural forms despite the molecular simplicity [53]. A recent ab-initio study by Brinkmann [57] describes formic acid dimerization as a stepwise process. The study found that the acyclic dimer, which is thermodynamically less stable, converts to the global minimum, a cyclic structure with two O–H···O–C hydrogen bonds. This cyclic dimer has a direct dissociation energy of 67.4 kJ mol⁻¹ and is 26.4 kJ mol⁻¹ more stable than the acyclic dimer [57]. Farfan et al. [54] used DFT to determine the equilibrium geometries of various formic acid dimer combinations (cis + cis, cis + trans, and trans + trans). These authors identified 21 distinct minima within a 69.0 kJ mol⁻¹ energy window on the potential energy surface for formic acid dimers.

In the gas phase, the dimer is more stable by ca. 58.6 kJ mol⁻¹ compared to the monomer [61]. The population of the cyclic H-bonded dimer configuration, dominant in the gas phase, is estimated to be around 23% [56,62], with approximately 95% of the vapor existing in the dimerized form at room temperature and normal pressure leading to non-ideal gas behavior [34]. Unlike most carboxylic acids, which retain dimeric structures in their crystalline form, formic acid exhibits an infinite “polymeric” crystal structure [63,64] – H-bonded chains with alternating O–H···O and C–H···O bonds [53,65,66]. Additionally, at high pressure, crystalline formic acid can adopt a cis molecular structure, differing from the trans structure found under normal conditions (where “trans” and “cis” refer to the positions of the hydrogen atoms) [53,67]. The liquid structure of formic acid is complex and consists of a hydrogen-bonded network with various chain-like structures, differing from both the gas and crystalline forms [53,68–71]. Small clusters with O–H···O bonds are held together by weaker C–H···O bonds, forming large branched structures, as suggested by recent first-principles and classical molecular dynamics simulations [53,62]. In addition to polymeric chains and cyclic dimers, an acyclic “open” dimer or a mixture of these species have been proposed as the main components of liquid formic acid [55]. The structure and dynamics of liquid formic acid related to H-bonding have been investigated using reverse Monte Carlo simulations by Jedlovszky et al. [69] Jedlovszky, Turi et al. performed Monte Carlo [58,68] and molecular dynamics (MD) [72] simulations using a potential model based on the ab-initio potential energy surface of formic acid dimers [58]. Hermida Ramón and Ríos [49] developed a polarizable force field for trans-formic acid for MD simulations in the liquid phase.

In this work, a multi-scale methodology, including force field-based Monte Carlo simulations, Quantum Mechanics, and equations of state is proposed to model the thermodynamic properties of a formic acid reacting system to explore its potential as a new reactive fluid in

thermodynamic cycles. This study is organized as follows: in Section 2, we provide technical details of the computational and mathematical methods. The methodologies for conducting two types of umbrella sampling Monte Carlo (MC) simulations—the dimer counter method and the potential of mean force method—are described, along with the thermodynamic modeling of the formic acid reactive mixture. In Section 3, we present and discuss the results. The thermodynamic properties are computed from Monte Carlo simulations and compared to Quantum Mechanics calculations. The dimer counter and PMF methods showed strong agreement, resulting in the enthalpies of HCOOH dimerization of ca. –60.46 kJ mol⁻¹ and ca. –62.91 kJ mol⁻¹, respectively. The QM-calculated enthalpy of dimerization is –60.48 kJ mol⁻¹, deviating by 0.02% from the dimer counter method and 4% from the PMF method. The QM and umbrella sampling MC simulations demonstrate high accuracy in reproducing the experimental thermodynamic properties of HCOOH dimerization. The Global Phase Equilibrium of the system, vaporization enthalpy, phase composition, vapor and liquid densities of the coexisting phases, and entropy as a function of temperature were obtained from the thermodynamic model and compared with Monte Carlo simulations in the Gibbs ensemble. The MC simulations agree well with the model, confirming the accuracy of the simulations in capturing the phase behavior of the system. Our findings are summarized in Section 4.

## 2. Methodology

### 2.1. Monte Carlo simulations

The dimerization of formic acid in thermodynamic equilibrium is simulated using the Monte Carlo (MC) Software Brick-CFCMC [73,74], an open-source molecular simulation code for the calculation of phase and reaction equilibria using state-of-the-art force field-based MC simulations in different ensembles, such as the NVT, NPT, grand-canonical, reaction, and the Gibbs ensemble [74]. To sample the dimerization reaction, it is crucial to explore all possible distances between the two molecules of formic acid in the Monte Carlo simulation in the canonical ensemble. Due to the strong interaction between the two monomers through hydrogen bonds, once a dimer forms, it tends to persist in each simulation cycle without dissociating. This phenomenon leads to a simulation dominated by dimers, resulting in an insufficient sampling of monomers at various distances, which hinders the accurate sampling of thermodynamic averages, such as the equilibrium constant and enthalpy of dimerization. To address this issue, we implement umbrella sampling to Brick-CFCMC, a technique that biases the sampling process to include a range of distances between the two molecules. By applying this method, we can obtain a more comprehensive representation of the system, enabling precise calculation of equilibrium properties. Additionally, while not employed in this study, the Aggregation-Volume-Bias Monte Carlo (AVBMC) moves [75, 76] have the potential to enhance configurational exploration by efficiently sampling system-wide volume changes. This approach could improve sampling efficiency in certain scenarios and would also be
worth exploring in future studies.

In a Monte Carlo simulation in the canonical ensemble, ensemble averages are computed by sampling configurations with a probability proportional to the Boltzmann factor ($e^{-\beta U(r^N)}$) as follows [77]:

$$ \langle A \rangle = \frac{\int dr^N A(r^N) \exp[-\beta U(r^N)]}{\int dr^N \exp[-\beta U(r^N)]} \qquad (2) $$

where $\beta = \frac{1}{k_B T}$, and $U(r^N)$ is the potential energy of a configuration. The conventional Metropolis acceptance rule for a trial move is given by [77,78]:

$$ \text{acc}(\text{old} \to \text{new}) = \min (1, \exp[-\beta \cdot (U_{\text{new}} - U_{\text{old}})]) \qquad (3) $$

where $U_{old}$ is the potential energy of the current configuration, and $U_{new}$ is the potential energy of the new configuration. In the umbrella sampling method, a simulation is performed in an ensemble $\pi$, in which a weight function $W(d(r^N))$ is added to the partition function, and the configurations are sampled with a probability proportional to [77]:

$$ \pi(r^N, d) = \exp[-\beta U(r^N) + W(d(r^N))] \qquad (4) $$

where $W(d(r^N))$ is a weight function (or biasing function) of the configuration that depends on the distance between the centers of two molecules $d$. The new acceptance rule for translation trial moves is given by:

$$ \text{acc}(\text{old} \to \text{new}) = \min (1, \exp[\mathcal{W}_{\text{new}}(d_{\text{new}}) - \mathcal{W}_{\text{old}}(d_{\text{old}}) - \beta \cdot (U_{\text{new}} - U_{\text{old}})]) \qquad (5) $$

where $d_{old}$ and $d_{new}$ are the old and new distances between dimers, respectively. The ensemble average $\langle A \rangle$ in the canonical ensemble can be computed by [79]:

$$ \langle A \rangle = \frac{\langle A \exp[-W] \rangle_{\pi}}{\langle \exp[-W] \rangle_{\pi}} \qquad (6) $$

where $\langle \dots \rangle_{\pi}$ is an ensemble average in the ensemble $\pi$ (Eq. (4)). This assumes that a single coordinate (distance) is sufficient to describe dimerization.

To sample all possible distances between the two molecules, a histogram of $W(d)$ must be obtained. This ensures that all distances are accommodated in simulation bins. Initially, $W(d)$ is set to zero. During a short simulation, the following iterative rule is applied:

$$ W(d) \to W(d) - \frac{1}{2} \ln(\text{Hist}(d)) \qquad (7) $$

where Hist(d) is the histogram count of observations for a given distance $d$ between the centers of two molecules. The procedure involves a series of short simulations, where the updated $W(d)$ from each simulation is used as the input for the next. The process continues until the resulting histogram of $W(d)$ becomes sufficiently flat, ensuring that all distances in the simulation box are adequately sampled. The final function $W(d)$ is then used in a longer simulation of the dimerization reaction. From this simulation, the equilibrium constants of dimerization can be obtained using two routes: the dimer counter method and the potential of mean force (PMF) method. Both routes are implemented in the umbrella sampling simulations, as the same set of sampled configurations produced in a single MC simulation can be used to compute the equilibrium constant of dimerization using either the dimer counter or the PMF method.

The dimer counter method uses geometric criteria from a study by Turner et al. [52] to identify whether the configuration consists of one dimer or two monomers of formic acid:

$$ 0.20 \text{ nm} \le d_{\text{Ofa2...Ofa1}} \le 0.34 \text{ nm} \qquad (8) $$

$$ 0.14 \text{ nm} \le d_{\text{Ofa2...Hfa2}} \le 0.24 \text{ nm} \qquad (9) $$

Both conditions must be met simultaneously to identify a dimer. The schematic representation of the geometric criteria for HCOOH dimer with the atoms labeled is visualized in Fig. 2.

The number of monomers ($N_{monomer}$) and dimers ($N_{dimer}$) are computed using the probabilities of observing two monomers ($p$) or one dimer ($1-p$):

$$ N_{\text{monomer}} = 2p \qquad (10) $$

$$ N_{\text{dimer}} = 1 - p \qquad (11) $$

The equilibrium constant $K$ is defined as:

$$ K = \frac{\frac{P_{\text{dimer}}}{P_0}}{\left(\frac{N_{\text{monomer}} \cdot RT}{P_0 \cdot V}\right)^2} = \frac{\frac{N_{\text{dimer}} \cdot RT}{P_0 \cdot V}}{\left(\frac{N_{\text{monomer}} \cdot RT}{P_0 \cdot V}\right)^2} = \frac{N_{\text{dimer}}}{N_{\text{monomer}}^2} \cdot \frac{P_0 \cdot V}{RT} \qquad (12) $$

where $P_0$ is the reference pressure (1 bar), $V$ is the volume of the simulation box, $R$ is the universal gas constant, and $T$ is the temperature in the studied range from 250 K to 300 K. Substituting $N_{monomer}$ and $N_{dimer}$ from Eqs. (10) and (11) into Eq. (12), the equilibrium constant is given as:

$$ K = \frac{(1-p)}{(2p)^2} \cdot \frac{P_0 \cdot V}{RT} \qquad (13) $$

The enthalpy of dimerization is determined using the Van't Hoff equation [80]:

$$ \frac{d \ln K}{d\left(\frac{1}{T}\right)} = -\frac{\Delta d H}{R} \qquad (14) $$

The second approach for computing the equilibrium constants of HCOOH dimerization involves the potential of mean force (PMF). Umbrella sampling is used for computing the PMF, which allows to understand the energy changes involved in the association/dissociation processes of noncovalent ligand-receptor pairs [81]. This method is particularly useful because it can be applied to a wide range of systems, from small molecular systems [82,83] to complex biological systems [81,84,85]. By computing the free energy along the dissociation pathway, the PMF provides thermodynamic details for molecular recognition. The PMF is obtained from umbrella sampling simulations using:

$$ PMF = - \ln \left( \frac{\text{Hist}(d) \cdot \exp(-W(d))}{\text{Norm}_2 \cdot V_{\text{bin}}} \right) \qquad (15) $$

where Hist(d) is the histogram count of observations for a given distance bin $d$, $W(d)$ is the weight function for the distance bin $d$, Norm$_2$ is a normalization factor for Hist(d), and $V_{bin}$ is the volume of the bin in the distance histogram, accounting for the spacing between discrete bins.

The change in Helmholtz free energy $\Delta F$ is calculated by integrating the PMF over the distance $d$, from an initial distance $d_0$ of ca. 2.6 Å to ca. 4.6 Å using:

$$ \Delta F = RT \cdot \int_{d_0}^{d} PMF(d')dd' \qquad (16) $$

This distance range is selected based on the probability distribution $p(d)$ of the distance between the centers of two HCOOH molecules computed using the geometrical criteria. Substituting the enthalpy and Helmholtz free energy change from Eqs. (17) and (18) into Eq. (19), the change in Gibbs free energy $\Delta G$ was determined:

$$ \Delta H = \Delta U + \Delta(PV) \qquad (17) $$

$$ \Delta F = \Delta U - T\Delta S \quad (\text{at } T = \text{const.}) \qquad (18) $$

$$ \begin{aligned} \Delta G &= \Delta H - T\Delta S = \Delta U + \Delta(PV) - T\Delta S \\ &= \Delta F + \Delta(PV) = \Delta F + RT\Delta n \quad (\text{at } T, V = \text{const.}) \end{aligned} \qquad (19) $$

where $\Delta n = 1$ (representing the dissociation of a dimer). The equilibrium constant $K$ is obtained from Eq. (20), and is used to calculate the enthalpy of dimerization from the Van't Hoff equation (Eq. (14)).

$$ K = \exp \left[ -\frac{\Delta G}{RT} \right] \qquad (20) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Graphical representation of the geometric criteria for the HCOOH dimer [52] with the atoms labeled.</figcaption></figure>

The molar entropy for HCOOH dimerization was calculated at 298.15 K for both the dimer counter method and the PMF method by:

$$ \Delta_d S = \frac{\Delta_d H - \Delta G}{298.15 \text{ K}} \qquad (21) $$

The composition of the vapor and liquid phases at saturation conditions is computed using the Gibbs Ensemble method at constant total volume, combined with the CFC method [73]. This ensemble consists of two simulation boxes that can exchange molecules. In the Gibbs ensemble, thermodynamic equilibrium of dimerization is achieved without the need for umbrella sampling, as both the formation of dimers and their subsequent dissociation into monomers are observed due to molecule exchanges between simulation boxes. Initially, each simulation box contains 200 molecules, and both boxes are identical in terms of the number of molecules and volume. Phase equilibrium densities are simulated for a temperature range of 335 K to 560 K. The mole fractions of monomers in each phase are computed using the implemented monomer counter based on the inverse dimer geometric criterion (where pairs of molecules that do not meet the dimer criterion defined in Eqs. (8) and (9), are considered monomers). For the temperature range from 335 K to 460 K, additional NPT simulations of the gas phase were performed using densities obtained from Gibbs Ensemble simulations. This approach was used to determine the mole fractions of HCOOH monomers in the gas phase with higher accuracy. Saturated vapor pressures are estimated by assuming an ideal gas phase as [86]:

$$ P_{\text{sat}} = \rho_{\text{liq}} \cdot k_{\text{B}} T \cdot \exp \left( \frac{\mu_{\text{liq}}^{\text{ex}}}{k_{\text{B}} T} \right) \qquad (22) $$

where $\rho_{\text{liq}}$ and $\mu_{\text{liq}}^{\text{ex}}$ are the number density and the excess chemical potential of the liquid phase from Gibbs Ensemble simulations, respectively. Despite the assumption of an ideal gas phase, this approach provides better agreement with experimental data compared to the computation from the series of NPT simulations of vapor phase, as shown in a study of Wasik et al. [45] The P-(x, y) and T-(x, y) diagrams are obtained and compared with results from the Equation of State. The enthalpy of vaporization can be calculated using the following equation for a non-ideal gas [87]:

$$ \Delta_{\text{vap}} H = (U_{\text{vap}} + P_{\text{sat}} V_{\text{vap}}) - (U_{\text{liq}} + P_{\text{sat}} V_{\text{liq}}) \qquad (23) $$

where $U_{\text{vap}}$ and $U_{\text{liq}}$ represent the potential energies of the vapor and liquid phases, respectively, R is the gas constant, $P_{\text{sat}}$ is the saturation pressure, $V_{\text{vap}}$ and $V_{\text{liq}}$ represent the molar volumes of the vapor and liquid phases, respectively. The entropy of the liquid and vapor phases

was calculated from [88]:

$$ s_i = \frac{U_i - \mu N_i + k_B T \ln P_{\text{sat}}}{T} \qquad (24) $$

where $U_i$ is the potential energy of the vapor or liquid phase, $\mu_i$ is the chemical potential of the vapor or liquid phase, $N_i$ is the number of molecules in the vapor or liquid phase, $k_B$ is the Boltzmann constant, T is the temperature, and $P_{\text{sat}}$ is the pressure at saturation. To check the correlation between the vapor pressure and the enthalpy of vaporization prediction, the Clausius–Clapeyron equation was used [89]:

$$ \Delta_{\text{vap}} H = \frac{dP}{dT} \cdot T \cdot V_{\text{vap}} = \frac{s_{\text{vap}} - s_{\text{liq}}}{V_{\text{vap}}} \cdot T \cdot V_{\text{vap}} \qquad (25) $$

The dimerization of HCOOH is studied using the ‘FF-2’ variant of the OPLS/AA force field for HCOOH [90], previously applied in our earlier research [45,47,48]. This force field is applied to the HCOOH molecule model, which was constructed and optimized at the B3LYP/6-31G(d) level of theory in a study of Wasik et al. [45] The HCOOH model is rigid, with point charges assigned to all atoms. The LJ parameters and partial charges for HCOOH used in this work are listed in Table S1 of the Supporting Information. The geometry of the HCOOH molecule is presented in Table S2 of the Supporting Information. The HCOOH force field has been validated by successful reproduction of the vapor–liquid equilibrium coexistence curve, saturated vapor pressures, and densities at various temperatures [45]. Additionally, the force field was used in our previous studies on the adsorption of HCOOH in M-MOF-74 [48], as well as the CO₂ hydrogenation reaction to HCOOH in UiO-66, Cu-BTC, IRMOF-1 [47], and M-MOF-74 [48], where it was applied together with the non-polarizable CO₂ and H₂ force fields [91]. The presence of dimers was confirmed in the gas phase of pure HCOOH and HCOOH/H₂O/NaCl simulations using the HCOOH force field, leading to non-ideal gas behavior and an inaccurate description of the HCOOH/H₂O azeotrope [45].

In all MC simulations, intermolecular interactions are modeled using Coulombic and Lennard-Jones (LJ) potentials. The Lorentz–Berthelot mixing rules [92] are applied for interactions between different LJ sites. The cutoff radius for intermolecular interactions is set to 10 Å, with interactions truncated and analytic tail corrections applied. Periodic boundary conditions are applied in all three directions. Electrostatic interactions in terms of energy are computed using the Ewald summation method [93], with the number of k-vectors in each direction ($K_{\text{max}}$) set to 8, and the damping parameter ($\alpha$) set to 0.32 Å⁻¹. The Ewald summation parameters correspond to a relative precision of 10⁻⁶.

The weight functions $W(d)$ for HCOOH dimerization equilibrium are obtained from a series of ten short consecutive umbrella sampling
simulations for each temperature within the range of 250 K to 300 K. These simulations involve 2 molecules of HCOOH in a simulation box with a length of 80 Å. Each simulation comprises 2·10² initialization MC cycles, 2·10⁵ equilibration MC cycles, and 10⁶ production MC cycles. A single MC cycle consists of N MC trial moves, where N is the total number of molecules at the start of the simulation. The final weight functions are used as input for extended umbrella sampling simulations of 2 molecules of HCOOH in a simulation box of the same size and at the corresponding temperatures to compute equilibrium constants. These simulations involve 2·10³ initialization MC cycles, 10⁶ equilibration MC cycles, and 10⁸ production MC cycles. The probabilities for selecting trial moves were 50% for translations and 50% for rotations. The Gibbs Ensemble simulations are performed with 2·10³ initialization MC cycles, 2·10⁵ equilibration MC cycles, and 10⁶ production MC cycles. The probabilities for selecting trial moves in these simulations were 28.57% for translations, 28.57% for rotations, 14.29% for λ changes, and 28.57% for CFC hybrid trial moves, which combine swap and identity changes [73]. To calculate the standard deviations of the computed values, all sets of simulations were conducted five times, starting from independent configurations and using different random number seeds. Random initial configurations were generated using the algorithms provided in Brick-CFCMC [73,74].

## 2.2. Thermodynamic modeling of formic acid formation

One of the successful tentative approaches used in the modeling of the thermodynamic properties of carboxylic acids is based on the representation of these systems as reactive mixtures [94,95]. It is possible to model the breaking (dissociation) and formation (association) of hydrogen bonds between, respectively, dimers and monomers, as a chemical reaction effect (see Fig. 1).

The theoretical model which has been applied in this work to characterize the formic acid system is similar to the one that has recently been used to represent the reactive N₂O₄ ⇌ 2NO₂ mixture [13], where the use of a cubic equation of state has already proven its successful application. This model is summarized below. For more details, the reader is referred to the description by Lasala et al. [13].

The equation of state which has been selected is the cubic Peng-Robinson equation of state (Eq. (26)), coupled with advanced mixing rules, EoS+a<sub>res</sub><sup>E,y</sup> (Eq. (27)) [96].

$$P(T, v, z) = \frac{RT}{v - b_m} - \frac{a_m}{v(v + b_m) + b_m(v - b_m)} \quad (26)$$

where T is the temperature, R is the universal gas constant, v is the molar volume, z is the vector of the molar composition, and a<sub>m</sub> and b<sub>m</sub> are the mixing energy and covolume parameters calculated as follows:

$$\left\{ \begin{array}{l} b_m = \sum_{i=1}^{NC} z_i b_i \\ \\ \frac{a_m}{b_m} = \sum_{i=1}^{NC} z_i \frac{a_i}{b_i} + \frac{a_{res}^{E,y}}{\Lambda_{EoS}} \end{array} \right. \quad (27)$$

where EoS+a<sub>res</sub><sup>E,y</sup> is the residual part of an excess Helmholtz energy model calculated from an activity coefficient (γ) model, Λ<sub>EoS</sub> is a numerical parameter being a function of the considered equation of state, the pure component energy and co-volume parameters, i.e. a<sub>i</sub> and b<sub>i</sub>, respectively, are calculated using the standard PR-78 equation of state [97] as a function of the critical temperature, critical pressure and acentric factor of the two components of the mixture (the monomer and the dimer of formic acid). The composition z is a result of the chemical equilibrium condition, the zero-Gibbs energy condition:

$$\Delta_R G(T, P, z) = 0 \quad (28)$$

In this work, the athermal version of Eq. (27) has been considered, i.e. the residual excess Helmholtz energy is null (a<sub>res</sub><sup>E,y</sup> = 0). This choice is motivated by the fact that the strong physical interactions related to hydrogen bonds are modeled here as chemical interactions instead of physical association effects.

To calculate the thermodynamic-phase equilibrium and energetic properties of the reactive system formed by monomers and dimers of, for example, formic acid with a cubic equation of state, it is necessary to determine some basic properties for the single molecules forming the mixtures: the monomeric and the dimeric form of, in this case, formic acid. These basic properties are the ideal gas thermochemical properties (standard enthalpy of formation, standard absolute entropy) and some thermophysical ones (the ideal gas heat capacity, critical temperature, critical pressure, and acentric factor). To calculate these properties, we have performed Quantum Mechanics (QM) calculations to determine the ideal gas properties, and then an optimization step aimed to fit the critical coordinates of the two pure components in order to match the experimental vaporization enthalpy and saturation pressures of the reactive mixture. Since the reactive system is at chemical equilibrium, it is not possible to experimentally determine the critical parameters of both the monomer and the dimer.

Quantum mechanics calculations were performed at the CBS-QB3 level of theory, with the Gaussian09 software [98], to obtain the lowest energy conformer. A post-processing of the results has been performed with GPOP [99] to derive thermodynamic properties. The results of these calculations are shown in where the expression for the ideal gas isobaric specific heat capacity is given by:

$$c_{P,l}^{ig}(T) = c_0 + c_1 \left( \frac{c_2/T}{\sinh(c_2/T)} \right)^2 + c_3 \left( \frac{c_4/T}{\cosh(c_4/T)} \right)^2 \quad (29)$$

We have then optimized the critical temperature, T<sub>c1</sub> and T<sub>c2</sub>, and the critical pressure, P<sub>c1</sub> and P<sub>c2</sub>, of the monomer (subscript 1) and of the dimer (subscript 2) by minimizing the following objective function:

$$\min_{\substack{T_{c1}, T_{c2} \\ P_{c1}, P_{c2}}} F = \frac{1}{n_{\Delta_{\text{vap}} H^{\text{exp}}}} \sum_{i=1}^{n_{\Delta_{\text{vap}} H^{\text{exp}}}} \left( \frac{\Delta_{\text{vap}} H^{\text{calc}} - \Delta_{\text{vap}} H^{\text{exp}}}{\Delta_{\text{vap}} H^{\text{exp}}} \right)^2 + \frac{1}{n_{P_{\text{VLE}}^{\text{exp}}}} \sum_{i=1}^{n_{P_{\text{VLE}}^{\text{exp}}}} \left( \frac{P_{\text{VLE}}^{\text{calc}} - P_{\text{VLE}}^{\text{exp}}}{P_{\text{calc}}^{\text{exp}}} \right)^2 + \left( \frac{T_c^{\text{calc}} - T_c^{\text{exp}}}{T_c^{\text{exp}}} \right)^2 + \left( \frac{P_c^{\text{calc}} - P_c^{\text{exp}}}{P_c^{\text{exp}}} \right)^2 \quad (30)$$

The experimental data used for this optimization are the ones reported for formic acid by the DIPPR database [100] (79 points for P<sub>VLE</sub> and 23 points for Δ<sub>vap</sub>H, see Tables S3 and S4 of the Supporting Information, respectively) and the ones obtained in this study. It is worth highlighting that the available data are quite dispersed and cover, for the vapor-liquid equilibrium pressure, a quite limited temperature range (up to 400 K). At a first but acceptably accurate approach, the acentric factor of dimers can be considered to be double the acentric factor of monomers. In this study, the acentric factor is assumed to be 0.1 for the monomer of formic acid and 0.2 for its dimeric form, based on the values proposed for ethane and butane [101], respectively.

## 3. Results and discussion

To investigate the dimerization of HCOOH at thermodynamic equilibrium, we first computed the weight functions W(d) and the observed probability distributions p(d) (normalized Hist(d)) of the distance between the centers of two HCOOH molecules from the umbrella sampling MC simulations at 250–300 K. Fig. 3 shows examples of W(d) and p(d) computed at 250 K. In Fig. 3a and b, the initial weight function (resulting from the first short consecutive umbrella sampling simulation) and the corresponding probability distribution of the distance between the centers of two HCOOH molecules are shown, respectively. In the initial W(d), the distances longer than ca. 4.6 Å are not sampled as the corresponding probability distribution p(d) is close to 0. The distances from ca. 2.9 Å to 4.6 Å, corresponding to dimer formation, are sampled with a very high probability. Since monomers are not sampled in the simulation, the equilibrium of HCOOH dimerization cannot be reproduced. In Fig. 3c and d, the final weight function (resulting from the last short consecutive umbrella sampling simulation) and the corresponding
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Weight function $W(d)$ in units of $k_B T$ and probability distribution $p(d)$ of the distance between the centers of two HCOOH molecules simulated at 250 K: (a) the weight function resulting from the first of a series of ten short umbrella sampling MC simulations, (b) the probability distribution of the observed distance between two molecules resulting from the first of a series of ten short umbrella sampling MC simulations, (c) the weight function resulting from the tenth of a series of ten short umbrella sampling MC simulations, (d) the probability distribution of the distance between two molecules resulting from the tenth of a series of ten short umbrella sampling MC simulations. The probabilities for the distance distribution sum up to 100%.</figcaption></figure>

observed probability distribution of the distance between the centers of two HCOOH molecules are shown, respectively. The distances are sampled starting from ca. 2.6 Å, which is the shortest distance between two HCOOH molecules possible due to intermolecular interactions. The distance range between two HCOOH molecules considered a dimer (2.6–4.6 Å) is sampled with a slightly higher probability than longer distances assigned to monomers (longer than 4.6 Å) by ca. 0.06. Due to the small difference between the probabilities of sampling dimers and monomers, the final weight function $W(d)$ was selected as the input for the extended umbrella sampling simulation of two HCOOH molecules at the corresponding temperature to determine the thermodynamic properties of HCOOH dimerization. The procedure used to compute $W(d)$ at 250 K was applied consistently for simulations performed at 260, 270, 290, and 300 K.

The equilibrium constants of HCOOH dimerization were computed from both routes of umbrella sampling MC simulations: dimer counter and potential of mean force. To obtain the equilibrium constants from the PMF method, it is necessary to calculate the change in Helmholtz free energy and subsequently the Gibbs free energy, see Eq. (20). The Helmholtz free energy change $\Delta F$ is calculated using Eq. (16) by integrating the PMF over the distance $d$, from an initial distance $d_0$ of ca. 2.6 Å to ca. 4.6 Å. In this range, two molecules are considered dimers according to the probability distribution $p(d)$ of the distance between the centers of two HCOOH molecules computed using the geometrical criteria [52] (Eqs. (8) and (9)). This approach ensures compatibility between the dimer counter and PMF methods and maintains consistency in the criteria used to identify a dimer. The example of PMF as a function of the distance between the centers of two HCOOH molecules is shown in Fig. 4.

The logarithms of the equilibrium constants for HCOOH dimerization computed using the dimer counter method and the potential of mean force (PMF) method are compared in Fig. 5. The changes in Gibbs free energy and equilibrium constants used to obtain the Arrhenius plots are listed in Table S5 of the Supporting Information. The values of $\ln K$ from both methods show close agreement, with an average deviation of ca. 2%. The logarithms of the equilibrium constants of HCOOH dimerization computed from MC methods differ from the experimental study of Buettner and Maurer [102] by approximately 2 units, while showing similar slopes of the linear trendlines. Despite the discrepancy in the absolute magnitude of the equilibrium constants potentially impacting the species distribution, the model accurately captures the qualitative behavior of the dimerization process. This shows the ability of the model to represent non-ideal interactions and trends in the system, which offers valuable insights for applications such as the design of heat pumps and other thermodynamic cycles. Compared to literature data from the study by Chao and Zwolinski [60] based on the statistical thermodynamic method, the values of $\ln K$ from the dimer counter and PMF methods are shifted by approximately 2 units, which slightly affects the slopes of the linear trendlines and consequently the enthalpy of dimerization. Compared to the computational study of formic acid dimerization in a carbon dioxide solvent by Turner et al. [52], the values of $\ln K$ from both methods are shifted by approximately 4 units. The enthalpies of HCOOH dimerization, derived from the slopes of the linear trendlines in the Arrhenius plots ($-\Delta_d H/R$), are provided in Table 1. The $\Delta_d H$ values differ by only 4% between the dimer counter and PMF methods, resulting in -60.46 kJ mol⁻¹ and -62.91 kJ mol⁻¹, respectively. The difference between the molar entropy of dimerization between the two routes is slightly higher at 7%. The deviation between the enthalpy of dimerization from Monte
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Potential of mean force (PMF) in units of $k_B T$ from an umbrella sampling MC simulation of two HCOOH molecules at 250 K as a function of the distance between the centers of two HCOOH molecules *d*. The value of Helmholtz free energy change $\Delta F$ calculated for this example by integrating the PMF over the distance *d*, from an initial distance $d_0$ of ca. 2.6 Å to ca. 4.6 Å, is -28.33 kJ mol$^{-1}$.</figcaption></figure>

Carlo simulations and the study by Chao and Zwolinski [60] is only ca. 3%. The deviation between the enthalpy of dimerization from Monte Carlo simulations and the computational study by Turner et al. [52] is more significant, at an average of 23%. In the study of formic acid dimerization by Turner et al. [52], carbon dioxide is present in the simulated system as a solvent, and a different HCOOH force field is used (from the study of Jedlovszky and Turi [58]). Due to these differences, umbrella sampling MC simulations show better agreement with experimental data of pure formic acid than the computations of Turner et al. [52] The dimerization enthalpies obtained from the dimer counter and PMF methods deviate by ca. 4% from P-V-T methods [60], ca. 9% from spectroscopic methods [60], and ca. 4% from the enthalpy calculated using experimental equilibrium constants [102]. The dimer counter method shows a smaller deviation than the PMF method, with an average difference of 4% from the experimental methods, while PMF results in an 8% difference. In sharp contrast, $\Delta_d H$ from the study of Turner et al. [52] deviates by ca. 15% from P-V-T methods, ca. 11% from spectroscopic methods, and ca. 15% from the enthalpy calculated using experimental equilibrium constants [102]. These results indicate that the umbrella sampling MC simulations used in our work provide an accurate reproduction of the enthalpy of HCOOH dimerization. The dimerization entropies obtained from MC methods also show good agreement with the value calculated using experimental equilibrium constants [102] deviating by only ca. 5%. Additionally, the results from the MC simulations were compared to the association model for dimerization of alkanoic acids from the study of Nagata et al. [103]. The logarithms of the equilibrium constants for the association of HCOOH molecules calculated from the association model differ by approximately 3 units from both MC results, leading to a significant deviation of about 31% in the enthalpy of dimerization (-42.4 kJ mol$^{-1}$), and the entropy of dimerization approximately three times smaller. The enthalpy of dimerization calculated from the association model does not reproduce the experimental methods accurately showing a deviation of ca. 27%. While association models are widely used and effective for many systems, these models rely on simplifying assumptions regarding molecular interactions, which can limit their accuracy in describing thermodynamic properties. In contrast, the methodology applied in this work provides a detailed molecular-level understanding, enabling the direct computation of dimerization equilibrium constants, enthalpy, and entropy without relying on pre-assumed association parameters.

The thermodynamical properties computed from MC simulations were also compared to the Quantum Mechanics calculations. The results of the Quantum Mechanics calculations are reported in. The enthalpy of monomer formation at 298.15 K calculated from QM (-379.07 kJ mol$^{-1}$) is in excellent agreement with the literature data reported by Chao and Zwolinski [60] (-378.57 kJ mol$^{-1}$), resulting in 0.1% difference. Similarly, the calculated enthalpy of dimer formation at 298.15 K (-818.617 kJ mol$^{-1}$) deviates by only 0.3% from the literature value of -820.94 kJ mol$^{-1}$ [60]. The entropies of monomer and dimer formation at 298.15 K obtained from QM are 254 J mol$^{-1}$K$^{-1}$ and 370 J mol$^{-1}$K$^{-1}$, respectively. This results in 2% of deviation from the reported entropy of monomer formation (248.88 J mol$^{-1}$K$^{-1}$) [60] and 11% from the reported entropy of dimer formation (332.67 J mol$^{-1}$K$^{-1}$) [60]. The logarithms of the equilibrium constants for HCOOH dimerization from QM calculations were compared to the dimer counter method and the PMF method in Fig. 5. Compared to MC methods, the ln K values from QM are shifted by approximately 0.16 units and result in the enthalpy of dimerization equal to -60.48 kJ mol$^{-1}$, deviating from the dimer counter method by 0.02% and from PMF by 4% ( Table 1). The value of $\Delta_d H$ calculated from QM differs similarly from the experimental data [60,102] as the MC simulations, resulting in an average deviation of ca. 4% from experimental methods, and the deviation of ca. 5% from the value calculated using the enthalpies of monomer and dimer formation at 298.15 K from the study by Chao and Zwolinski [60]. The molar entropy of dimerization calculated from QM (-138 J mol$^{-1}$K$^{-1}$) differs by an average of 3% from both MC methods and is in good agreement with the value calculated using experimental equilibrium constants (-150.29 J mol$^{-1}$K$^{-1}$) [102] deviating by ca. 8%, as well as with the value obtained from the entropies of monomer and dimer formation at 298.15 K (-165.09 J mol$^{-1}$K$^{-1}$) [60] deviating by ca. 16%.

The QM and umbrella sampling MC simulations used in this study have demonstrated a high degree of accuracy in reproducing the experimental thermodynamic properties of HCOOH dimerization. Our results show that the dimer counter and PMF methods are in excellent agreement. Notably, the dimer counter method provides a more precise match with experimental data compared to the PMF method.

The results of the Quantum Mechanics calculations () and the optimization procedure (Table 3) have enabled plotting the Global Phase Equilibrium Diagram of the system (Fig. 6), the vaporization enthalpy
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Arrhenius plot for HCOOH dimerization at 250–300 K computed from Quantum Mechanics and both routes of umbrella sampling MC simulations: dimer counter and potential of mean force. The logarithms of the equilibrium constants are compared with literature data from the experimental study by Buettner and Maurer [102], the study by Chao and Zwolinski [60] using a statistical thermodynamic method, the computational study by Turner et al. [52], and the association model by Nagata et al. [103]. The enthalpies of dimerization are calculated from slopes of the linear trendlines. The error bars are smaller than the size of the symbols.</figcaption></figure>

**Table 1**

Enthalpy and entropy of dimerization computed from Quantum Mechanics and both routes of umbrella sampling MC simulations: dimer counter and potential of mean force. The subscripts show uncertainties computed as the standard deviation from five independent simulations. The computed results are compared with literature data [52,60,102,103].

<table><thead><tr><th>Method</th><th>Δ<sub>d</sub>H/<br>(kJ/mol)</th><th>Δ<sub>d</sub>S/<br>(J/mol/K)</th></tr></thead><tbody><tr><td>MC - Dimer counter</td><td>-60.46<sub>0.45</sub></td><td>-137.36<sub>1.76</sub></td></tr><tr><td>MC - PMF</td><td>-62.91<sub>0.64</sub></td><td>-146.98<sub>2.27</sub></td></tr><tr><td>QM</td><td>-60.48</td><td>-138</td></tr><tr><td>Buettner and Maurer [102]</td><td>-59.29</td><td>-150.29</td></tr><tr><td>Chao and Zwolinski [60]</td><td>-63.81</td><td>-165.09</td></tr><tr><td>Turner et al. [52]</td><td>-50.3</td><td></td></tr><tr><td>Nagata et al. [103]</td><td>-42.44</td><td></td></tr><tr><td>P-V-T [60]</td><td>-59.12</td><td></td></tr><tr><td>Spectroscopy [60]</td><td>-56.57</td><td>-46.99</td></tr></tbody></table>

**Table 2**

Results of QM calculations for the monomer and dimer of HCOOH, including the enthalpy and entropy of monomer and dimer formation at 298.15 K, as well as the expression for the ideal gas isobaric specific heat capacity (Eq. (29)).

<table><thead><tr><th>Compound</th><th>Δ<sub>f</sub>H°<sub>i</sub><sup>-298.15K</sup>/<br>(kJ/mol)</th><th>Δ<sub>f</sub>S°<sub>i</sub><sup>-298.15K</sup>/<br>(J/mol/K)</th><th>c<sub>0</sub>/<br>(J/kmol/K)</th><th>c<sub>1</sub>/<br>(J/kmol/K)</th></tr></thead><tbody><tr><td>HCOOH</td><td>-379.07</td><td>254</td><td>36 322.4</td><td>46 443.9</td></tr><tr><td>(HCOOH)<sub>2</sub></td><td>-818.617</td><td>370</td><td>77 376.5</td><td>97 286.3</td></tr></tbody></table><table><thead><tr><th>Compound</th><th>c<sub>2</sub>/<br>K</th><th>c<sub>3</sub>/<br>(J/kmol/K)</th><th>c<sub>4</sub>/<br>K</th></tr></thead><tbody><tr><td>HCOOH</td><td>724.28</td><td>41 856.3</td><td>2168.01</td></tr><tr><td>(HCOOH)<sub>2</sub></td><td>712.16</td><td>92 502.6</td><td>1979.58</td></tr></tbody></table>

(Fig. 7), and the non-isothermal but unique phase diagram showing the molar composition of the liquid phase (x₁) and of the vapor phase (y₁) as a function of the pressure (and thus of the temperature) of this monovalent system (Fig. 8). Additionally, the vapor and liquid densities of the coexisting phases of the reactive HCOOH mixture, as a function of temperature (Fig. 9) was calculated, as well as T-s diagram (Fig. 10). More details on the methodology are provided in Lasala et al. [13]. The results of Monte Carlo simulations in the Gibbs ensemble are shown in

**Table 3**

Results of the optimization procedure for the monomer and the dimer of HCOOH.

<table><thead><tr><th>Compound</th><th>T<sub>c</sub>/K</th><th>P<sub>c</sub>/bar</th><th>Acentric Factor</th></tr></thead><tbody><tr><td>HCOOH</td><td>350</td><td>40</td><td>0.1</td></tr><tr><td>(HCOOH)<sub>2</sub></td><td>633</td><td>38</td><td>0.2</td></tr></tbody></table>

Table S6 of the Supporting Information.

Fig. 6 shows the Global Phase Equilibrium Diagram that highlights the complex behavior of the monomer and dimer forms of formic acid, as well as the critical phenomena that arise from the interactions. The results of thermodynamic model calculations are compared to MC simulations in the Gibbs ensemble and the experimental data [100]. The figure features the saturation pressure of the monomer of formic acid, the saturation pressure of the dimer of formic acid, the reactive mixture vapor-liquid equilibrium pressure curve, and the locus of the critical points of the inert mixtures formed by monomers and dimers. As temperature increases from 281.5 K (which is the melting point of formic acid [104]) to a critical value of 350 K, the saturation pressure of the formic acid monomer also rises, reflecting the typical behavior of a pure substance where higher temperatures require higher pressures to maintain the substance in a liquid state. The curve starts at a low-pressure value of ca. 10 bar and gradually increases with temperature to a critical pressure of ca. 40 bar. The curve of the saturation pressure of the HCOOH dimer lies below the monomer saturation pressure curve, indicating that the dimer has a lower vapor pressure compared to the monomer at the same temperature. This is consistent with the fact that dimers, being larger and having stronger intermolecular forces, are less volatile than monomers. The curve also shows an upward trend as temperature increases from 281.5 K to a critical value of 633 K, where the critical pressure is as high as 38 bar. The vapor-liquid equilibrium pressure of the reactive mixture includes contributions from both monomers and dimers. This curve lies between the monomer and dimer saturation pressure curves, reflecting the combined effects of both species in the equilibrium. As temperature rises from 281.5 K to ca. 570 K, the equilibrium pressure also increases up to ca. 65 bar, determined by the ongoing reaction between monomers and dimers.
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Global Phase Equilibrium Diagram, showing the saturation pressure of the monomer of formic acid (orange curve), the saturation pressure of the dimer of formic acid (yellow curve), the reactive mixture vapor–liquid equilibrium pressure curve (blue curve), and the locus of the critical points of the inert mixtures formed by monomers and dimers (violet curve), calculated with the thermodynamic model: (a) in the temperature range from 281.5 K to 633 K, and (b) a zoom-in on the low-temperature range. The green points represent experimental vapor–liquid equilibrium pressures [100] while the blue cross represents the experimental critical coordinate of the equilibrium mixture. The black points represent vapor–liquid equilibrium pressures computed from MC simulations in the Gibbs ensemble. Vapor–liquid equilibrium properties of HCOOH computed from MC simulations in the Gibbs ensemble are listed in Table S6 of the Supporting Information.</figcaption></figure>

The results are compared with vapor–liquid equilibrium pressures computed using MC simulations and experimental data [100], showing reasonable agreement across the examined conditions and confirming the thermodynamic model prediction. The locus of the critical points for the inert mixture of monomers and dimers shows a parabolic shape, peaking at ca. 490 K and ca. 80 bar, and then declining. This behavior is characteristic of critical phenomena for systems of type I according to the classification of Van Konynenburg and Scott [105–107].

The vaporization enthalpy calculated from the thermodynamic model is compared to MC simulations in the Gibbs ensemble and the experimental data [100] in Fig. 7. In the temperature region from 281.5 K

to 500 K, the vaporization enthalpy is relatively high (ca. 450 kJ mol⁻¹), as indicated by both the experimental data and the thermodynamic model. The curve remains almost flat or slightly increases with temperature. This suggests that the amount of heat required to vaporize the substance does not change significantly with temperature in this range. The MC simulations show a significant deviation in this region, predicting higher values of vaporization enthalpy than observed experimentally and by the model. This discrepancy may be due to an overestimation of the hydrogen bonding strength resulting from the HCOOH force field. Overestimating hydrogen bonding leads to a more stabilized liquid phase, potentially resulting in a too-high
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Vaporization enthalpy as a function of temperature of the reactive mixture at chemical equilibrium. The blue line represents calculations from the thermodynamic model, while the green points represent the experimental data [100]. The black points represent the results of MC simulations in the Gibbs ensemble.</figcaption></figure>

concentration of molecular clusters, in the simulation. The force field used for formic acid has not been analyzed for these higher-order clusters present in the liquid and vapor phases, as the mole fraction was only calculated for monomers. If the simulations predict an excessive number of trimers or higher-order clusters compared to reality, this would increase the cohesive energy of the liquid, leading to a higher vaporization enthalpy. Similar issues are presented in the study of Schnabel et al. [59], using simulation results of Chialvo et al. [108], where a large deviation of the computed enthalpy of vaporization from experimental data is shown in the low-temperature region. Additionally, the NVT simulations of two HCOOH molecules showed that without implementing the umbrella sampling technique, there is a related issue: due to the strong interactions between monomers via hydrogen bonds, once a dimer forms, it tends to persist in each simulation cycle without dissociating. This persistence further suggests that the hydrogen bond strength may be overestimated for the force field. In the high-temperature region, particularly above 500 K, the curve shows a sharp decline in the enthalpy of vaporization. Above 500 K, formic acid becomes unstable, which is reflected in the decreasing enthalpy of vaporization observed in experiments. This indicates that the heat required to vaporize the mixture becomes significantly less as the temperature approaches its critical point, where the liquid and vapor phases become indistinguishable. This behavior is typical as the thermal energy of the system increases, making it easier for molecules to escape the liquid phase. The experimental data and the thermodynamic model start to converge with the MC simulations in this temperature range, suggesting a better agreement between theoretical and experimental observations. All three data sets – experimental, thermodynamic model, and MC simulations – show this sharp decrease, which is consistent with the behavior expected near the critical point. The better agreement at higher temperatures might indicate that the effects of overestimated hydrogen bonding and cluster formation are reduced as thermal fluctuations dominate and the system approaches the critical point. The enthalpy of vaporization calculated from the MC simulation data using Eq. (23) was checked for consistency with the Clausius–Clapeyron equation [89]. The comparison resulted in a good agreement between the two calculation methods with a deviation of only ca. 6%, showing the reliability of the applied methodology.

**Fig. 8a** shows a non-isothermal phase diagram for a reactive formic acid mixture at chemical equilibrium, with the pressure *P* as a function of the mole fractions *x*₁ and *y*₁, representing the HCOOH monomers in the liquid and vapor phases, respectively. The phase behavior calculated from the thermodynamic model is compared with the results from MC simulations in the Gibbs ensemble. In the liquid phase, the mole fraction of the monomers increases rapidly from *x*₁ = 0 to ca. 0.35 with pressure up to ca. 65 bar, which suggests that more monomers are retained in the liquid phase, reducing the tendency to vaporize. The results of MC simulations align closely with the thermodynamic model curve, showing the accurate reproduction of the behavior of the liquid phase. This agreement indicates that the force field used in the simulations is well-suited for modeling the liquid-phase behavior of HCOOH monomers. The vapor-phase curve behaves differently, initially increasing more slowly at low pressures and low mole fractions. It eventually reaches a maximum mole fraction of ca. 0.5, where the curve bends and begins to decrease. This behavior is characteristic of a vapor phase where, beyond a certain pressure, further increases lead to the condensation of monomers, thereby reducing the mole fraction in the vapor. The MC simulation data for the vapor phase also match the thermodynamic model curve, particularly in the higher pressure range of 20–40 bar. At pressures lower than 20 bar, there is a slight deviation between the simulation data points and the model curve, but the overall trend is consistent. This slight discrepancy at lower pressures could be due to the limitations in the MC simulations, such as the challenge of accurately capturing the complex phase behavior in reactive mixtures at chemical equilibrium. The phase diagram shows that the thermodynamic model is successful in describing the non-isothermal phase behavior of the HCOOH mixture, with the liquid and vapor phase curves exhibiting expected behaviors as pressure varies. The MC simulations agree well with the model across most of the pressure and mole fraction range, confirming the validity of the simulations in capturing the phase behavior of the system. The phase compositions are also presented as a function of temperature in Fig. 8b, where the trends in the mole fraction of HCOOH monomers exhibit similar patterns.

**Fig. 9** shows vapor and liquid densities of the coexisting phases of the reactive HCOOH mixture, as a function of temperature calculated
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Phase diagrams of the reactive HCOOH mixture at chemical equilibrium: (a) non-isothermal *P*-{*x*, *y*}, and (b) *T*-{*x*, *y*}. The curves represent the mole fraction of HCOOH monomer in the liquid phase (blue curve) and in the vapor phase (red curve), calculated with the thermodynamic model. The points represent the mole fraction of HCOOH monomer in the liquid phase (blue points) and the vapor phase (red points), resulting from MC simulations in the Gibbs ensemble.</figcaption></figure>

from the thermodynamic model, compared to MC simulation in the Gibbs ensemble. The vapor phase density, starting from very low densities and increasing as the temperature approaches the critical point, is in good agreement with MC simulations and closely matches the

experimental data, especially at lower densities. However, the liquid phase density of the thermodynamic model is underestimated by ca. 400 kg m⁻³ at 335 K but still follows the general trend. As temperature increases, the liquid phase curve calculated from the thermodynamic
<figure><img src="image_12.png" /><figcaption>Fig. 9. Vapor and liquid densities of the coexisting phases of the reactive HCOOH mixture at chemical equilibrium, as a function of temperature. The curves represent the liquid phase density (blue curve) and the vapor phase density (red curve), calculated with the thermodynamic model. The points represent the liquid phase density (blue points) and the vapor phase density (red points), resulting from MC simulations in the Gibbs ensemble. The green curve represents experimental data [100].</figcaption></figure>

model slowly converges to MC simulation and experimental data. While the MC simulations agree well with experimental data, effectively capturing the experimental liquid phase behavior up to the critical point, there is still potential for improving the thermodynamic model for formic acid, particularly in accurately reproducing the liquid phase density. The low accuracy of cubic equations of state in reproducing the liquid phase density is well known. Nevertheless, an improvement in modeling liquid densities could be achieved with the inclusion of a volume-translation term in the cubic equation of state [109].

Fig. 10 shows T-s diagram of the reactive HCOOH mixture at chemical equilibrium, calculated with the thermodynamic model and compared to MC simulations in the Gibbs ensemble. The liquid phase has a relatively low entropy compared to the vapor phase at the same temperature, with the highest difference at 281.5 K (ca. 2.9 kJ mol⁻¹K⁻¹ and ca. 4.5 kJ mol⁻¹K⁻¹ for liquid and vapor phases, respectively). This difference is due to the structured nature of the liquid phase compared to the disordered vapor phase. In the vapor phase, temperature increases significantly from 281.5 K to 500 K while entropy changes relatively little (from ca. 4.5 kJ mol⁻¹K⁻¹ to 4.7 kJ mol⁻¹K⁻¹). In this phase, additional heat primarily increases the temperature, rather than further increasing entropy due to a low heat capacity. In sharp contrast, in the liquid phase temperature increases more steadily from 281.5 K to ca. 580 K with an increase in entropy from ca. 2.9 kJ mol⁻¹K⁻¹ to ca. 4.5 kJ mol⁻¹K⁻¹. This gradual change reflects a higher heat capacity of the liquid phase. The results of MC simulations in the Gibbs ensemble are in good agreement with the thermodynamic model, reproducing the trend in entropy of the liquid and vapor phases. The values of entropy of the liquid and vapor phase are shifted from the thermodynamic model by ca. 1.5 kJ mol⁻¹K⁻¹ and ca. 1 kJ mol⁻¹K⁻¹, respectively. The close alignment between the thermodynamic model and MC simulations confirms the accuracy of the multi-scale methodology in capturing the behavior of both phases for a wide range of temperatures.

## 4. Conclusions

This work introduces a multi-scale methodology that integrates force field-based Monte Carlo simulations, Quantum Mechanics, and equations of state to characterize the thermodynamics of the formic

acid dimerization reaction. This approach explores the potential application of formic acid as a new reactive fluid in thermodynamic cycles, particularly in the design of heat pumps. Accurate knowledge of thermodynamic properties is crucial for optimizing heat transfer and energy efficiency in these systems. To compute thermodynamic properties, such as the equilibrium constants, enthalpy, and entropy of dimerization, umbrella sampling was implemented in Monte Carlo simulations. This technique biases the sampling process to include a range of distances between the two molecules. The logarithms of the equilibrium constants for HCOOH dimerization were computed at 250–300 K using both the dimer counter and the potential of mean force methods. Both methods showed strong agreement, with an average deviation of only 1.8%. The enthalpies of HCOOH dimerization derived from the Arrhenius plots were –60.46 kJ mol⁻¹ and –62.91 kJ mol⁻¹ for the dimer counter and PMF methods, respectively. The molar entropy of dimerization computed from the dimer counter method is –137.36 J mol⁻¹K⁻¹, while from the PMF method, it is –146.98 J mol⁻¹K⁻¹. The results of Monte Carlo simulations were compared to Quantum Mechanics calculations. The QM-calculated enthalpy of dimerization is –60.48 kJ mol⁻¹, deviating by 0.02% from the dimer counter method and 4% from the PMF method. The QM-calculated molar entropy of dimerization is –138 J mol⁻¹K⁻¹, differing by an average of 3% from both MC methods. The QM and umbrella sampling MC simulations employed in this study demonstrated high accuracy in reproducing the experimental thermodynamic properties of HCOOH dimerization, with the dimer counter method providing the most precise match with experimental data. To calculate the thermodynamic phase equilibrium properties, the cubic Peng–Robinson equation of state, coupled with an athermal version of advanced mixing rules, was applied. Ideal gas properties, which cannot be experimentally measured, were determined using Quantum Mechanics calculations. These properties are necessary for an optimization step aimed at fitting the critical coordinates of the two pure components to match the experimental vaporization enthalpy and saturation pressures of the reactive mixture. The Global Phase Equilibrium of the system, vaporization enthalpy, phase composition, vapor and liquid densities of the coexisting phases as a function of temperature, and entropy as a function of temperature were obtained from the thermodynamic model and compared with Monte Carlo simulations
<figure><img src="image_13.png" /><figcaption>Fig. 10. T-s diagram of the reactive HCOOH mixture at chemical equilibrium, calculated with the thermodynamic model. The blue curve represents the liquid phase, where temperature increases gradually with entropy. The red curve represents the vapor phase, where temperature increases significantly while entropy changes relatively little due to a low heat capacity. The points represent the entropy of the liquid phase (blue points) and the vapor phase (red points), resulting from MC simulations in the Gibbs ensemble.</figcaption></figure>

in the Gibbs ensemble. Overall, the MC simulations agree well with the model, especially for saturation pressure and the non-isothermal phase diagram of the HCOOH mixture, confirming the validity of the simulations in capturing the phase behavior of the system. These results are directly relevant to the design of heat pumps using formic acid as a working fluid, where precise control over thermodynamic properties is essential for optimizing performance. However, it is important to note that this study does not consider the kinetics of the dimerization reaction, which could influence the performance of heat pumps. Future work could incorporate reaction kinetics into the simulations, potentially through kinetic Monte Carlo or Molecular Dynamics simulations, to provide a more comprehensive understanding of the system behavior. To conclude, the presented methodology has proven effective in accurately determining the chemical equilibrium properties of 2HCOOH ⇌ (HCOOH)₂. It could be adapted to preliminarily predict the thermodynamic properties of other similar reactive systems, such as other carboxylic acids, thereby broadening its applicability in the design of heat pumps and other thermodynamic devices. It is important to note that our rigid HCOOH force field was not optimized or fitted to the dimerization reaction. As a first step, the Monte Carlo simulation only considers monomers and dimers. A more realistic description would include trimers and more complex structures.

### CRediT authorship contribution statement

**Dominika O. Wasik:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. **Silvia Lasala:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Olivier Herbinet:** Writing – review & editing, Software, Methodology, Formal analysis, Data curation. **Konstantin Samukov:** Methodology, Investigation. **Sofía Calero:** Supervision, Resources, Funding acquisition. **Thijs J.H. Vlugt:** Writing – review & editing, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This research has been supported by the Eindhoven Institute for Renewable Energy Systems (EIRES). S.L. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon Europe research and innovation program (grant agreement No. 101040994).

### Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.fluid.2025.114356.

### Data availability

Data will be made available on request.

### References

[1] European Commission, Cluster 5: Climate, energy and mobility, 2023, https://ec.europa.eu/info/research-and-innovation/funding/funding-opportunities/funding-programmes-and-open-calls/horizon-europe/cluster-5-climate-energy-and-mobility_en. (Access date: 5 June 2024).

[2] IEA, CO₂ emissions in 2023, Paris, 2024, https://www.iea.org/reports/CO₂-emissions-in-2023; (Accessed date: 15 May 2024).

[3] A. Barakat, S. Lasala, P. Arpentinier, J.-N. Jaubert, The original and impactful exploitation of chemical energy in heat pumps, Chem. Eng. J. Adv. 12 (2022) 100400.

[4] G. Kosmadakis, Estimating the potential of industrial (high-temperature) heat pumps for exploiting waste heat in EU industries, Appl. Therm. Eng. 156 (2019) 287-298.

[5] D. van de Bor, C.I. Ferreira, A.A. Kiss, Low grade waste heat recovery using heat pumps and power cycles, Energy 89 (2015) 864-873.
[6] A. Khanlari, A. Sözen, B. Sahin, G.D. Nicola, F. Afshari, Experimental investigation on using building shower drain water as a heat source for heat pump systems, Energy Sources A: Recover. Util. Environ. Eff. (2020) 1–13.

[7] IEA, Heat pumps, paris, 2021, https://www.iea.org/energy-system/buildings/heat-pumps. (Access date: 5 May 2024).

[8] C. Forman, I.K. Muritala, R. Pardemann, B. Meyer, Estimating the global waste heat potential, Renew. Sustain. Energy Rev. 57 (2016) 1568–1579.

[9] K. Zhang, Y. Shen, C. Duwig, Identification of heat transfer intensification mechanism by reversible $N_2O_4$ decomposition using direct numerical simulation, Int. J. Heat Mass Transfer 182 (2022) 121946.

[10] C. Li, H. Wang, Power cycles for waste heat recovery from medium to high temperature flue gas sources – from a view of thermodynamic optimization, Appl. Energy 180 (2016) 707–721.

[11] Z. Xu, R. Wang, C. Yang, Perspectives for low-temperature waste heat recovery, Energy 176 (2019) 1037–1043.

[12] S. Lasala, R. Privat, O. Herbinet, P. Arpentinier, D. Bonalumi, J.-N. Jaubert, Thermo-chemical engines: Unexploited high-potential energy converters, Energy Convers. Manage. 229 (2021) 113685.

[13] S. Lasala, K. Samukov, H.M. Polat, V. Lachet, O. Herbinet, R. Privat, J.-N. Jaubert, O.A. Moultos, K. De Ras, T.J.H. Vlugt, Application of thermodynamics at different scales to describe the behaviour of fast reacting binary mixtures in vapour-liquid equilibrium, Chem. Eng. J. 483 (2024) 148961.

[14] A. Barakat, S. Lasala, P. Arpentinier, P. Tobaly, J.-N. Jaubert, Understanding the thermodynamic effects of chemically reactive working fluids in the stirring engine, Energy Convers. Manag.: X 22 (2024) 100573.

[15] A. Krasin, V. Nesterenko, Dissociating gases: a new class of coolants and working substances for large power plants, At. Energy Rev. 9 (1971) 177–194.

[16] H.M. Huang, R. Govind, Use of dissociating gases in brayton cycle space power systems, Ind. Eng. Chem. Res. 27 (1988) 803–810.

[17] A. Sorokin, Dissociating nitrogen tetroxide ($N_2O_4$) as a working fluid in thermodynamic cycles, Nucl. Sci. Eng. 72 (1979) 330–346.

[18] R.J. Stochl, Potential Performance Improvement Using a Reacting Gas (Nitrogin Tetroxide) as the Working Fluid in a Closed Brayton Cycle, 1979, p. 1060.

[19] M. Binotti, C.M. Invernizzi, P. Iora, G. Manzolini, Dinitrogen tetroxide and carbon dioxide mixtures as working fluids in solar tower plants, Sol. Energy 181 (2019) 203–213.

[20] A. Belkadi, F. Lovell, V. Gerbaud, L. Vega, Modeling the vapor–liquid equilibrium and association of nitrogen dioxide/dinitrogen tetroxide and its mixtures with carbon dioxide, Fluid Phase Equilib. 266 (2008) 154–163.

[21] L.E.S. de Souza, U.K. Deiters, Modeling of the $N_2O_4$–$NO_2$ reacting system, Phys. Chem. Chem. Phys. 2 (2000) 5606–5613.

[22] G. Manzolini, M. Binotti, D. Bonalumi, C. Invernizzi, P. Iora, $CO_2$ mixtures as innovative working fluid in power cycles applied to solar plants, Techno-Econ. Assess. Sol. Energy 181 (2019) 530–544.

[23] D.N. Seshadri, D.S. Viswanath, N.R. Kuloor, Thermodynamic properties of the system $N_2O_4 \rightleftharpoons 2NO_2 \rightleftharpoons 2NO + O_2$, AIChE J. 16 (1970) 420–425.

[24] D.F. Stai, F. Bizjak, S.E. Stephanou, Thermodynamic properties of nitrogen tetroxide, J. Spacecr. Rockets 2 (1965) 742–745.

[25] E. Bourrasseau, V. Lachet, N. Desbiens, J.-B. Maillet, J.-M. Teuler, P. Ungerer, Thermodynamic behavior of the $CO_2 + NO_2/N_2O_4$ mixture: A Monte Carlo simulation study, J. Phys. Chem. B 112 (2008) 15783–15792.

[26] S. Lasala, Reactive fluids for intensified thermal energy conversion, Proj. Repos. J. 13 (2022) 102–105, https://www.europeandissemination.eu/article/reactive-fluids-for-intensified-thermal-energy-conversion/18808. (Access date: 05 May 2024).

[27] S. Lasala, ERC-REACHER, 2022, https://www.univ-lorraine.fr/erc-reacher/. (Access date: 5 May 2024).

[28] R.I. Slavchov, J.K. Novev, S. Mosbach, M. Kraft, Vapor pressure and heat of vaporization of molecules that associate in the gas phase, Ind. Eng. Chem. Res. 57 (2018) 5722–5731.

[29] J. Chocholoušová, J. Vacek, P. Hobza, Potential energy and free energy surfaces of the formic acid dimer: Correlated ab initio calculations and molecular dynamics simulations, Phys. Chem. Chem. Phys. 4 (2002) 2119–2122.

[30] J. Eppinger, K.-W. Huang, Formic acid as a hydrogen energy carrier, ACS Energy Lett. 2 (2016) 188.

[31] S. Moret, P. Dyson, G. Laurenczy, Direct synthesis of formic acid from carbon dioxide by hydrogenation in acidic media, Nat. Commun. 5 (2014) 4017.

[32] Fact.MR, Formic acid market analysis by formic acid of 85% concentration, 90% concentration, 94% concentration, and 99% concentration from 2023 to 2033, 2023, https://www.factmr.com/report/4279/formic-acid-market; (Accessed date: 15 May 2024).

[33] M. Aresta (Ed.), Carbon Dioxide Recovery and Utilization, fourth ed., Dordrecht, Kluwer Academic Publishers, 2010.

[34] J. Hietala, A. Vuori, P. Johnsson, I. Pollari, W. Reutemann, H. Kieczka, Ullmann's Encyclopedia of Industrial Chemistry. Formic Acid, seventh ed., John Wiley & Sons, Ltd, 2016, 1–22.

[35] Fact.MR, Formic acid market analysis by concentration (formic acid of 85% concentration, 90% concentration, 94% concentration, 99% concentration), by application (formic acid for animal feed, agricultural products, leather tanning), by region - global insights 2021–2031, 2021, https://www.factmr.com/report/4279/formic-acid-market. (Access date: 20 August 2022).

[36] M. Rumayor, A. Dominguez-Ramos, A. Irabien, Formic acid manufacture: Carbon dioxide utilization alternatives, Appl. Sci. 8 (2018) 914.

[37] B.J. O'Neill, E.I. Gürbüz, J.A. Dumesic, Reaction kinetics studies of the conversions of formic acid and butyl formate over carbon-supported palladium in the liquid phase, J. Catalysis 290 (2012) 193.

[38] H.W. Gibson, Chemistry of formic acid and its simple derivatives, Chem. Rev. 69 (1969) 673.

[39] H. Yang, J.J. Kaczur, S.D. Sajjad, R.I. Masel, Electrochemical conversion of $CO_2$ to formic acid utilizing sustaininon™ membranes, J. CO₂ Util. 20 (2017) 208–217.

[40] M. Scibioh, B. Viswanathan, Electrochemical reduction of carbon dioxide: A status report, Proc. Indian Nat. Sci. Acad. 70 (2004) 407.

[41] X. Lu, D.Y.C. Leung, H. Wang, M.K.H. Leung, J. Xuan, Electrochemical reduction of carbon dioxide to formic acid, ChemElectroChem 1 (2014) 836–849.

[42] B. Hu, C. Guild, S.L. Suib, Thermal, electrochemical, and photochemical conversion of $CO_2$ to fuels and value-added products, J. CO₂ Util. 1 (2013) 18–27.

[43] M. Ramdin, A.R.T. Morrison, M. de Groen, R. van Haperen, R. de Kler, E. Irtem, A.T. Laitinen, L.J.P. van den Broeke, T. Breugelmans, J.P.M. Trusler, W.d. Jong, T.J.H. Vlugt, High-pressure electrochemical reduction of $CO_2$ to formic acid/formate: Effect of pH on the downstream separation process and economics, Ind. Eng. Chem. Res. 58 (2019) 22718–22740.

[44] M. Ramdin, A.R.T. Morrison, M. de Groen, R. van Haperen, R. de Kler, L.J.P. van den Broeke, J.P.M. Trusler, W. de Jong, T.J.H. Vlugt, High pressure electrochemical reduction of $CO_2$ to formic acid/formate: A comparison between bipolar membranes and cation exchange membranes, Ind. Eng. Chem. Res. 58 (2019) 1834–1847.

[45] D.O. Wasik, H.M. Polat, M. Ramdin, O.A. Moultos, S. Calero, T.J.H. Vlugt, Solubility of $CO_2$ in aqueous formic acid solutions and the effect of NaCl addition: A molecular simulation study, J. Phys. Chem. C 126 (2022) 19424–19434.

[46] M. Ghara, P.K. Chattaraj, A ca computational study on hydrogenation of $CO_2$, catalyzed by a bridged B/N frustrated Lewis pair, Struct. Chem. 30 (2019) 1067–1077.

[47] D.O. Wasik, A. Martín-Calvo, J.J. Gutiérrez-Sevillano, D. Dubbeldam, T.J.H. Vlugt, S. Calero, Enhancement of formic acid production from carbon dioxide hydrogenation using metal-organic frameworks: Monte Carlo simulation study, Chem. Eng. J. 467 (2023) 143432.

[48] D.O. Wasik, J.M. Vicent-Luna, S. Rezaie, A. Luna-Triguero, T.J.H. Vlugt, S. Calero, The impact of metal centers in the M-MOF-74 series on formic acid production, ACS Appl. Mater. & Interfaces 16 (2024) 45006–45019.

[49] J.M. Hermida Ramón, M.A. Rílós, A new intermolecular polarizable potential for cis-formic acid. Introduction of many-body interactions in condensed phases, Chem. Phys. 250 (1999) 155–169.

[50] W.L. Chameides, D.D. Davis, Aqueous-phase source of formic acid in clouds, Nature 304 (1983) 427–429.

[51] B. Wang, H. Hou, Y. Gu, New mechanism for the catalyzed thermal decomposition of formic acid, J. Phys. Chem. A 104 (2000) 10526–10528.

[52] C.H. Turner, Monte Carlo simulation of formic acid dimerization in a carbon dioxide solvent, J. Phys. Chem. B 108 (2004) 11716–11721.

[53] R. Chelli, R. Righini, S. Califano, Structure of liquid formic acid investigated by first principle and classical molecular dynamics simulations, J. Phys. Chem. B 109 (2005) 17006–17013.

[54] P. Farfán, A. Echeverri, E. Diaz, J.D. Tapia, S. Gómez, A. Restrepo, Dimers of formic acid: Structures, stability, and double proton transfer, J. Chem. Phys. 147 (2017) 044312.

[55] M. Gantenberg, M. Halupka, W. Sander, Dimerization of formic acid-an example of a noncovalent reaction mechanism, Chem. – A Eur. J. 6 (2000) 1865–1869.

[56] P. Rodziewicz, N.L. Doltsinis, Formic acid dimerization: Evidence for species diversity from first principles simulations, J. Phys. Chem. A 113 (2009) 6266–6274.

[57] N.R. Brinkmann, G.S. Tschumper, G. Yan, H.F. Schaefer, An alternative mechanism for the dimerization of formic acid, J. Phys. Chem. A 107 (2003) 10208–10216.

[58] P. Jedlovszky, L. Turi, A new five-site pair potential for formic acid in liquid simulations, J. Phys. Chem. A 101 (1997) 2662–2665.

[59] T. Schnabel, M. Cortada, J. Vrabec, S. Lago, H. Hasse, Molecular model for formic acid adjusted to vapor–liquid equilibria, Chem. Phys. Lett. 435 (2007) 268–272.

[60] J. Chao, B.J. Zwolinski, Ideal gas thermodynamic properties of methanoic and ethanoic acids, J. Phys. Chem. Ref. Data 7 (1978) 363–377.

[61] M.D. Taylor, J. Bruton, The vapor phase dissociation of some carboxylic acids. II. Formic propionic acids¹,², J. Am. Chem. Soc. 74 (1952) 4151–4152.

[62] I. Bakó, J. Hutter, G. Pálinkás, Car-Parrinello molecular dynamics simulation of liquid formic acid, J. Phys. Chem. A 110 (2006) 2188–2194.

[63] I. Nahringbauer, A Reinvestigation of the Structure of Formic Acid (at 98 K), Acta Crystallogr. B 34 (1978) 315–318.

[64] Z. Berkovitch-Yellin, L. Leiserowitz, Atom-atom potential analysis of the packing characteristics of carboxylic acids. A study based on experimental electron-density distributions, J. Am. Chem. Soc. 104 (1982) 4052–4064.
[65] B.L. Bhargava, Y. Yasaka, M.L. Klein, Hydrogen evolution from formic acid in an ionic liquid solvent: A mechanistic study by ab initio molecular dynamics, J. Phys. Chem. B 115 (2011) 14136-14140.

[66] L. Weng, C. Chen, J. Zuo, W. Li, Molecular dynamics study of effects of temperature and concentration on hydrogen-bond abilities of ethylene glycol and glycerol: Implications for cryopreservation, J. Phys. Chem. A 115 (2011) 4729-4737.

[67] D.R. Allan, S.J. Clark, Impeded dimer formation in the high-pressure crystal structure of formic acid, Phys. Rev. Lett. 82 (1999) 3464-3467.

[68] P. Jedlovszky, L. Turi, Role of the C-H···O hydrogen bonds in liquids: A Monte Carlo simulation study of liquid formic acid using a newly developed pair-potential, J. Phys. Chem. B 101 (1997) 5429.

[69] P. Jedlovszky, G.P.I. Bakó, J. Dore, Structural investigation of liquid formic acid, Mol. Phys. 86 (1995) 87-105.

[70] S. Nasr, M.-C. Bellissent-Funel, R. Cortès, X-ray and Neutron Scattering Studies of Liquid Formic Acid DCOOD at Various Temperatures and Under Pressure, J. Chem. Phys. 110 (1999) 10945-10952.

[71] L. Bellamy, R. Lake, R. Pace, Hydrogen bonding in carboxylic acids—II. Monocarboxylic acids, Spectrochim. Acta 19 (1963) 443-449.

[72] P. Mináry, P. Jedlovszky, M. Mezei, L. Turi, A comprehensive liquid simulation study of neat formic acid, J. Phys. Chem. B 104 (2000) 8287.

[73] R. Hens, A. Rahbari, S. Caro-Ortiz, N. Dawass, M.Erdős, A. Poursaeidesfahani, H.S. Salehi, A.T. Celebi, M. Ramdin, O.A. Moultos, D. Dubbeldam, T.J.H. Vlugt, Brick-CFCMC: Open source software for Monte Carlo simulations of phase and reaction equilibria using the continuous fractional component method, J. Chem. Inf. Model. 60 (2020) 2678-2682.

[74] H.M. Polat, H.S. Salehi, R. Hens, D.O. Wasik, A. Rahbari, F. de Meyer, C. Houriez, C. Coquelet, S. Calero, D. Dubbeldam, O.A. Moultos, T.J.H. Vlugt, New features of the open source Monte Carlo software brick-CFCMC: Thermodynamic integration and hybrid trial moves, J. Chem. Inf. Model. 61 (2021) 3752-3757.

[75] B. Chen, J.I. Siepmann, A novel Monte Carlo algorithm for simulating strongly associating fluids: Applications to water, hydrogen fluoride, and acetic acid, J. Phys. Chem. B 104 (2000) 8725-8734.

[76] B. Chen, J.I. Siepmann, Improving the efficiency of the aggregation-volume-bias Monte Carlo algorithm, J. Phys. Chem. B 105 (2001) 11275-11282.

[77] T.J.H. Vlugt, J.P.J.M. van der Eerden, M. Dijkstra, B. Smit, D. Frenkel, Introduction to Molecular Simulation and Statistical Thermodynamics, Delft, The Netherlands, Delft University of Technology, 2008.

[78] D. Frenkel, B. Smit, Understanding Molecular Simulation: from Algorithms to Applications, third ed., Elsevier, 2023.

[79] G. Torrie, J. Valleau, Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling, J. Comput. Phys. 23 (1977) 187-199.

[80] P. Atkins, J.De. Paula, J. Keeler, Atkins' Physical Chemistry, 12th ed., Oxford University Press, 2023.

[81] W. You, Z. Tang, C.-e.A. Chang, Potential mean force from umbrella sampling simulations: What can we learn and what is missed? J. Chem. Theory Comput. 15 (2019) 2433-2443, PMID: 30811931.

[82] N. Bansal, Z. Zheng, D.S. Cerutti, K.M. Merz, On the fly estimation of host-guest binding free energies using the movable type method: Participation in the SAMPL5 blind challenge, J. Comput. Aided Mol. Des. 31 (2017) 47-60.

[83] J. Yin, N.M. Henriksen, D.R. Slochower, M.R. Shirts, M.W. Chiu, D.L. Mobley, M.K. Gilson, Overview of the SAMPL5 host-guest challenge: Are we doing better? J. Comput. Aided Mol. Des. 31 (2017) 1-19.

[84] H. Sun, S. Tian, S. Zhou, Y. Li, D. Li, L. Xu, M. Shen, P. Pan, T. Hou, Revealing the favorable dissociation pathway of type II kinase inhibitors via enhanced sampling simulations and two-end-state calculations, Sci. Rep. 5 (2015) 8457.

[85] J.-K. Zhou, D.-Y. Yang, S.-Y. Sheu, The molecular mechanism of ligand unbinding from the human telomeric G-quadruplex by steered molecular dynamics and umbrella sampling simulations, Phys. Chem. Chem. Phys. 17 (2015) 12857-12869.

[86] N. Dawass, R.R. Wanderley, M. Ramdin, O.A. Moultos, H.K. Knuutila, T.J.H. Vlugt, Solubility of carbon dioxide, hydrogen sulfide, methane, and nitrogen in monoethylene glycol; Experiments and molecular simulation, J. Chem. Eng. Data 66 (2021) 524-534.

[87] J. Wang, H. Tingjun, Application of molecular dynamics simulations in molecular property prediction I: Density and heat of vaporization, J. Chem. Theory Comput. 7 (2011) 2151-2165.

[88] T.L. Hill, Thermodynamics of Small Systems, J. Chem. Phys. 36 (1962) 3182-3197.

[89] D. Koutsoyiannis, Clausius-Clapeyron equation and saturation vapour pressure: Simple theory reconciled with practice, Eur. J. Phys. 33 (2012) 295.

[90] F.J. Salas, E. Nunez-Rojas, J. Alejandro, Stability of formic acid/pyridine and isonicotinamide/formamide cocrystals by molecular dynamics simulations, Theor. Chem. Accounts 136 (2017) 1-12.

[91] D.O. Wasik, J.M. Vicent-Luna, A. Luna-Triguero, D. Dubbeldam, T.J.H. Vlugt, S. Calero, The impact of metal centers in the M-MOF-74 series on carbon dioxide and hydrogen separation, Sep. Purif. Technol. 339 (2024) 126539.

[92] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, second ed., Oxford University Press, 2017.

[93] B.A. Wells, A.L. Chaffee, Ewald summation for molecular simulations, J. Chem. Theory Comput. 11 (2015) 3684-3695.

[94] P. Grenzheuser, J. Gmehling, An equation of state for the description of phase equilibria and caloric quantities on the basis of the chemical theory, Fluid Phase Equilib. 25 (1986) 1-29.

[95] J. Gmehling, Ed, Chemical Thermodynamics for Process Simulation, second ed., Weinheim, Wiley-VCH, 2019.

[96] S. Lasala, P. Chiesa, R. Privat, J.-N. Jaubert, VLE properties of CO₂ – based binary systems containing N₂, O₂ and ar: experimental measurements and modelling results with advanced cubic equations of state, Fluid Phase Equilibria 428 (2016) 18-31.

[97] D.-Y. Peng, D. Robinson, The Characterization of the Heptanes and Heavier Fractions for the GPA Peng-Robinson Programs, Gas Processors Association, 1978.

[98] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheese-man, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, et al., Gaussian 09, Revision A1, Gaussian, Inc., Wallingford, CT, 2009.

[99] A. Miyoshi, GPOP software, 2002, http://akrmys.com/gpop/.

[100] DECHEMA, DIPPR 801 database, 2023, https://dechema.de/en/dippr801.html.

[101] C.A. Passut, R.P. Danner, Acentric factor. A valuable correlating parameter for the properties of hydrocarbons, Ind. Eng. Chem. Process. Des. Dev. 12 (1973) 365-368.

[102] R. Büttner, G. Maurer, Dimerisierung einiger organischer Säuren in der gasphase, Berichte der Bunsenges. Für Phys. Chem. 87 (1983) 877-882.

[103] I. Nagata, K. Gotoh, K. Tamura, Association model of fluids. Phase equilibria and excess enthalpies in acid mixtures, Fluid Phase Equilib. 124 (1996) 31-54.

[104] W.M. Haynes, CRC Handbook of Chemistry and Physics, 91st ed., Boca Raton, FL, CRC Press Inc., 2010.

[105] R.L. Scott, Models for phase equilibria in fluid mixtures, Acc. Chem. Res. 20 (1987) 97-107.

[106] P.H. van Konynenburg, Critical Lines and Phase Equilibria in Binary Mixtures (Ph.D. thesis), University of California, Los Angeles, USA, 1968.

[107] P.H. van Konynenburg, R.L. Scott, Critical lines and phase equilibriums in binary Van der Waals mixtures, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 298 (1980) 495-540.

[108] A.A. Chialvo, M. Kettler, I. Nezbeda, Effect of the range of interactions on the properties of fluids. 2. Structure and phase behavior of acetonitrile, hydrogen fluoride, and formic acid, J. Phys. Chem. B 109 (2005) 9736-9750.

[109] A. Peneloux, E. Rauzy, A consistent correction for Redlich-Kwong-Soave volumes, Fluid Phase Equilib. 8 (1982) 7-23.