# Computation of drug solvation free energy in supercritical CO2: alternatives
to all-atom computer simulations

Kalikin N.N  Budkov Y.A  Kolesnikov A.L  Ivlev D.V  Krestyaninov M.A  Kiselev
M.G  G.A. Krestov Institute of Solution Chemistry of the Russian Academy of
Sciences, Laboratory of NMR Spectroscopy and Numerical Investigations of
Liquids, Akademicheskaya Str. 1, 153045, Ivanovo, Russia  Tikhonov Moscow
Institute of Electronics and Mathematics, School of Applied Mathematics,
National Research University Higher School of Economics, 34, Tallinskaya
Ulitsa, 123458, Moscow, Russia  Institut für Nichtklassische Chemie e.V.,
Permoserstr. 15, 04318, Leipzig, Germany

###### Abstract

Despite the modern level of development of computational chemistry methods and
technological progress, fast and accurate determination of solvation free
energy remains a huge problem for physical chemists. In this paper, we
describe two computational schemes that can potentially solve this problem. We
consider systems of poorly soluble drug compounds in supercritical carbon
dioxide. Considering that the biggest contribution among all intermolecular
interactions is made by van der Waals interactions, we model solute and
solvent particles as coarse-grained ones interacting via the effective
Lennard-Jones potential. The first proposed approach is based on the classical
density functional theory and the second one relies on molecular dynamics
simulation of the Lennard-Jones fluid. Sacrificing the precision of the
molecular structure description while capturing the phase behavior of the
fluid with sufficient accuracy, we propose computationally advantageous paths
to obtaining the solvation free energy values with the accuracy satisfactory
for engineering applications. The agreement reached between the results of
such coarse-graining models and the experimental data indicates that the use
of the all-atom molecular dynamic simulations for the studied systems seems to
be excessive.

###### keywords:

Solvation free energy, statistical thermodynamics, classical density
functional theory, coarse-graining models, molecular dynamics simulation,
quantum chemistry, supercritical fluid.

## INTRODUCTION

Solvation free energy is an extremely important thermodynamic quantity for a
wide range of physico-chemical applications, including calculation of
solubility values, estimation of partition coefficients, infinite dilution
activity coefficients and other parameters. Nowadays there are a number of
computational methods developed to determine solvation free energy values.
Regarding the systems that include poorly soluble drugs and drug-like
compounds it is rather popular among scientists to rely on the results of the
all-atom molecular dynamics (AAMD) simulations based on alchemical free energy
methods [1, 2, 3, 4, 5, 6]. Such approaches model a set of nonphysical
intermediate states to compute the values of free energy of solute molecule
transfer from the gas phase to the solution. More details regarding different
alchemical methods can be found elsewhere [7]. Despite their evident
efficiency and similarly clear time expenditures and computational cost, these
methods can have other drawbacks. They include high sensitivity to the force
field choice and uncertainty arising from the choice of the partial atomic
charges computation method. Both of these aspects of the AAMD calculation can
lead to the inconsistency of the final solvation free energy values up to
several kcal/mol[8, 1, 9]. At the same time, the procedure of the force-field
parametrization is itself a nontrivial task and it is not always easy to
justify some of the approaches physically. Thus, if the aim is to determine
the free energy quickly, the utilization of the AAMD techniques is rather
questionable. On the other hand, as it has been recently shown [10, 11, 12,
13], the results of the coarse-graining models using united-atom force fields
are comparable in accuracy with those of the all-atom ones and in some cases
even outperform them. Another approach to the estimation of drug solvation
free energy in supercritical CO2 (scCO2), popular among chemical engineers, is
based on applying cubic equations of state, for instance, the Peng-Robinson
equation of state (PR EOS) for binary mixtures, which is a generalization of
the classical van der Waals equation of state [14, 15]. However, the PR EOS
contains several adjustable parameters, the determination of which requires an
experimental data set (solubility isotherms, for instance). That is why the PR
EOS is usually used to approximate the existing experimental data. Thus, in
order to apply the PR EOS for solvation free energy estimation, it is
necessary to perform preliminary experimental measurements that make this
approach much more complicated. The next step toward a more robust description
of the phase equilibrium is a family of SAFT (Statistical Associating Fluid
Theory) type models [16]. As they include the contributions of different
interaction terms and the parameters of these models hold a strict physical
meaning, they have been successfully utilized to correlate and predict the
solvation properties of a variety of solutes [17, 18, 19, 20, 21, 22].
Although there is a rather substantial decrease in the number of binary
interaction parameters, when multicomponent systems are considered, as
compared to the cubic EOS with advanced mixing rules, the whole
parametrization procedure is still cumbersome and the choice of the best one
is itself a nontrivial task [23]. High accuracy of the calculated solvation
free energy values with a error less than 1 kcal/mol for sufficiently large
data sets can be also achieved by using quantum mechanical calculations on
continuum solvents (COSMO-based models, SMD)[24, 25, 26, 27, 28]. The cost of
such accuracy, however, is the complex parametrization routine and the
necessity of its replication for new systems, if they are dramatically
different from the original learning set [29].

Quite promising is the method of solvation free energy estimation based on the
classical density functional theory (cDFT). In papers [30, 31], the authors
proposed a method for calculating hydration free energy of a set of
hydrophobic solutes. This method can be considered as a coarse-grained one,
since within its framework real molecules of both the solute and the solvent
are replaced by spherically symmetric particles interacting with each other
through the effective Lennard-Jones (LJ) potential. The utilization of the
Weeks-Chandler-Andersen (WCA) procedure for the LJ potential and application
of the fundamental measure theory [32] (FMT) to account for the short-range
hard core correlations allowed the authors to fit the pressure and the surface
tension of bulk water to the experimental values by changing the cutoff
radius. The latter, in turn, enabled them to describe hydration of a set of
hydrocarbons. Even though such cDFT approach remains a method of experimental
data approximation, it has distinct advantages over the methods based on the
implications of the PR EOS. The main one is that, in contrast to the PR EOS-
based approach [33], which in its nature is a variant of the mean-field theory
taking into account only long-range correlations of solution molecules, the
cDFT-based method accounts for the short-range correlations of the solvent
molecules due to the packing effects in the solute molecule solvation shell.
Moreover, within the cDFT approach, all the interaction parameters are related
directly to the effective LJ interaction potentials and, thereby, have a clear
physical meaning.

The development, in a sense, of the cDFT application for the solvation free
energy calculation can be traced by a number of papers [34, 35, 36, 37],
proposing the so-called molecular density functional theory approach, which
consists of several steps. The first one is computation of the site-site pair-
distribution functions from the AAMD simulation of the pure solvent. The
second one is obtaining of the direct correlation function from the Molecular
Ornstein–Zernike integral equation, and, finally, calculation of the energy
functional, the value of which, being the difference between the system grand
thermodynamic potential with and without the solute molecule, gives the
solvation free energy of the solute. Although such approach is faster than the
full-scale classical AAMD simulation, the first step is still dependent on the
force field choice. The results obtained for a series of halide ions are in
good agreement with the experiment, but the question is what results could be
achieved for drug-like molecules and how the calculation of such systems would
influence the computational speed.

In this paper, we discuss two coarse-graining techniques that can be used for
fast and sufficiently accurate estimation of solvation free energy values of
sparingly soluble drug compounds in a scCO2 medium. One of the potential
applications of the obtained data can be subsequent calculation of the
solubility of the compounds, which is a prerequisite for using supercritical
technologies, e.g. micronization or cocrystallization [38, 39], in order to
enhance the final aqueous solubility and thus the bioavailability of drugs.
The proposed approaches can make solvation free energy computation faster than
in the aforementioned techniques, preserving the accuracy of the obtained
values. The first one is based on cDFT, but, in contrast to the approaches
discussed above, the only input values required to compute the interaction
potential parameters are the critical parameters of the solvent and solute,
which can be often found in literature. The second approach is based on the
coarse-grained MD (CGMD) simulations of the LJ fluid, with the interaction
parameters determined on the basis on the law of the corresponding states,
knowing the critical parameters of the system compounds and critical
parameters of the LJ fluid.

## METHODOLOGY

## Classical Density Functional Theory

The details of the proposed density functional theory-based method are
presented in Appendix. Here we briefly outline the main points. Within such
approach CO2 molecules are coarse-grained to spherically symmetric particles
interacting through the effective pairwise LJ potential, the parameters of
which can be obtained by fitting the CO2 liquid-vapor critical point. Starting
from the expression of the grand thermodynamic potential for the CO2 fluid in
the external potential field $V_{ext}(\mathbf{r})$ created by the fixed solute
molecule located at the origin, one gets

$\Omega[\rho(\mathbf{r})]=F_{int}[\rho(\mathbf{r})]+\int
d\mathbf{r}\rho(\mathbf{r})(V_{ext}(\mathbf{r})-\mu),$ (1)

where $F_{int}[\rho(\mathbf{r})]$ is the intrinsic Helmholtz free energy of
the CO2 fluid and $\mu$ is the chemical potential of the bulk fluid at certain
temperature and pressure values. The intrinsic Helmholtz free energy, in its
turn, can be comprised of two contributions

$F_{int}[\rho(\mathbf{r})]=k_{B}T\int
d\mathbf{r}\rho(\mathbf{r})[\ln(\Lambda^{3}\rho(\mathbf{r}))-1]+F_{ex}[\rho(\mathbf{r})],$
(2)

where the first one is the Helmholtz free energy of the ideal gas and the
second one is the excess Helmholtz free energy of the fluid; $\Lambda$ is the
thermal de Broglie wavelength. One can write the total excess free energy in
the form

$F_{ex}[\rho(\mathbf{r})]=F_{hs}[\rho(\mathbf{r})]+F_{att}[\rho(\mathbf{r})].$
(3)

The hard spheres contribution (the first term on the right hand side) is
determined within Rosenfeld’s version of the fundamental measure theory (FMT)
[32] as follows

$F_{hs}[\rho(\mathbf{r})]=k_{B}T\int d\mathbf{r}\Phi(\mathbf{r}),$ (4)

where $\Phi(\mathbf{r})$ is the excess free energy density (expressed in
$k_{B}T$ units), which is the function of weighted densities. The attractive
contribution to the excess free energy is described within the mean-field
approximation [40]

$F_{att}[\rho(\mathbf{r})]=\frac{1}{2}\int d\mathbf{r}\rho(\mathbf{r})\int
d\mathbf{r}^{\prime}\rho(\mathbf{r}^{\prime})\phi_{WCA}(\mathbf{r}-\mathbf{r}^{\prime}),$
(5)

with the effective WCA pair potential of the attractive interactions.

The solute molecule is modeled as a particle creating an external LJ potential
with the effective parameters of the interaction between this molecule and the
molecules of the fluid, which are determined according to the Berthelot-Lorenz
mixing rules: $\sigma_{sf}=(\sigma_{ss}+\sigma_{ff})/2$ and
$\varepsilon_{sf}=\sqrt{\varepsilon_{ss}\varepsilon_{ff}}$. The total pressure
of the fluid bulk phase takes the following form

$P=\rho_{b}k_{B}T\frac{1+\eta+\eta^{2}}{(1-\eta)^{3}}+\frac{1}{2}B_{WCA}\rho_{b}^{2}.$
(6)

where $\eta=\pi d_{BH}^{3}\rho_{b}/6$ is the packing fraction of the hard
sphere system, the effective Barker-Henderson (BH) diameter is determined by
the Pade approximation
[41]:$d_{BH}=\sigma_{ff}(1.068\varepsilon_{ff}/k_{B}T+0.3837)/(\varepsilon_{ff}/k_{B}T+0.4293)$,
and the following auxiliary function corresponding to the attractive
contribution is introduced

$B_{WCA}=-\frac{32\sqrt{2}}{9}\pi\varepsilon_{ff}\sigma_{ff}^{3}+\frac{16}{3}\pi\varepsilon_{ff}\sigma_{ff}^{3}\left[\left(\frac{\sigma_{ff}}{r_{c}}\right)^{3}-\frac{1}{3}\left(\frac{\sigma_{ff}}{r_{c}}\right)^{9}\right].$
(7)

The parameters of the interaction potential between two molecules of CO2
($\sigma_{ff}$, $\varepsilon_{ff}$) can be obtained, as it was mentioned
above, by fitting the respective parameters of the liquid-gas critical point,
i.e. solving the following system of equations, using the known critical
parameters of CO2

$\frac{\partial P}{\partial\rho_{b}}\Bigr{|}_{\begin{subarray}{c}T=T_{c}\\\
\rho_{b}=\rho_{c}\end{subarray}}=0,~{}~{}\frac{\partial^{2}P}{\partial\rho_{b}^{2}}\Bigr{|}_{\begin{subarray}{c}T=T_{c}\\\
\rho_{b}=\rho_{c}\end{subarray}}=0.$ (8)

In the same manner one can find the parameters of the LJ potential for the
active compound ($\sigma_{ss}$, $\varepsilon_{ss}$). NIST [42] provides
accurate values of the critical parameters for CO2. At the same time, the
critical parameters of the solutes can be only estimated. The values of the
calculated potential parameters are presented in Table 1. The problems that
can potentially arise from the ambiguities of critical parameter estimation
are discussed in the next sections. After the iterative minimization of the
grand thermodynamic potential with respect to the density and finding its
equilibrium profile, one can obtain the solvation free energy as the excess
grand thermodynamic potential:

$\Delta G_{solv}=\Omega[\rho(\mathbf{r})]-\Omega[\rho_{b}].$ (9)

We would like to note that $\Omega[\rho_{b}]$ is calculated at
$V_{ext}(\mathbf{r})=0$.

Table 1: Parameters of the LJ potential of CO2 and drug compounds obtained from the fitting of the corresponding critical points using EOS implemented within the cDFT-based approach. compound |  $\sigma$, $[$\mathrm{\SIUnitSymbolAngstrom}$]$ |  $\varepsilon$, [K]
---|---|---
CO2 | 3.363 | 218.738
aspirin | 6.553 | 548.700
ibuprofen | 7.301 | 539.206
carbamazepine | 7.180 | 565.911

## All-atom MD simulation

For comprehensive investigation we have conducted a set of AAMD simulations
for the system of aspirin (ASP) in scCO2 to find out how the choice of the
AAMD simulation parameters affects the final solvation free energy values. We
have compared the results for three different force fields describing the ASP
molecule – GROMOS [43, 44], GAFF [45] and OPLS-AA [46], and two models
describing the CO2 molecules - Zhang’s [47] and TraPPE [48] ones. We have
established that the absolute values of the ASP solvation free energy obtained
using the GROMOS force field with the help of the Automated Topology Builder
system [43, 49, 50] are abnormally high when compared to those calculated by
the two other force-fields. As the procedure of accurate force field
parametrization for obtaining precise final values was not the goal of the
present study, rather, we wanted to test the performance of the generic force
fields without preliminary preparations, the results obtained with the GROMOS
force field will not be discussed.

The simulation details are similar to those we used in our previous studies
[51, 52]. Here we outline the main aspects. The AAMD simulations were
conducted in the Gromacs 4.6.7 program package [53, 54, 55, 56]
[http://www.gromacs.org]. The cell contained 1024 CO2 molecules and 1 ASP
molecule. The cell was constructed in the Packmol program [57], the CO2
density at different state parameters was taken from the NIST database [42].
The Bennett acceptance ratio (BAR) [58] approach was chosen to compute the
free energy surface. For each simulation we performed equilibration for 100 ps
with a 1 fs step and 500 ps with a 2 fs step in the NVT and NPT ensembles,
respectively, and the production run of 5 ns with the step of 2 fs.

We calculated partial atomic charges for the ASP molecule by the Merz-Kollman
method [59, 60] using Gaussian 09 software [61] with the PBE functional and
$\mathrm{6-311+}$ $\mathrm{+g(2d,p)}$ basis set. We averaged partial atomic
charges over two most stable ASP conformers (see Appendix). Such procedure
improved the agreement between the solvation free energy values of
carbamazepine (CBZ) obtained experimentally and computed within the AAMD
simulation in our previous paper [52].

For each computation of the solvation free energy surface we conducted 12
independent simulations, each corresponding to a different pair of coupling
parameters, describing the interactions between the solvent molecules and the
solute molecule. The potential function is linearly dependent on the coupling
parameters. The LJ solute-solvent interaction parameters were calculated from
the atomic parameters by the Berthelot-Lorentz mixing rules. We have chosen
the following set of alchemical coefficients $\\{\lambda_{LJ},\lambda_{C}\\}$:
$\\{0.0,0.0\\}$, $\\{0.2,0.0\\}$, $\\{0.5,0.0\\}$, $\\{1.0,0.0\\}$,
$\\{1.0,0.2\\}$, $\\{1.0,0.3\\}$, $\\{1.0,0.4\\}$, $\\{1.0,0.5\\}$,
$\\{1.0,0.6\\}$, $\\{1.0,0.7\\}$, $\\{1.0,0.8\\}$, $\\{1.0,1.0\\}$, where the
subscripts ”C” and ”LJ” correspond to the scaling parameters of the Coulomb
and Lennard-Jones potentials. The differences in the free energy between the
states with neighboring values of the coupling parameters were calculated
using the “g$\\_$bar” tool of the Gromacs package.

## Coarse-grained MD simulation

Besides the cDFT-based approach discussed above, we also propose a method of
total coarse-graining of the molecular structure of the sparingly soluble drug
compounds for accelerated computation of the solvation free energy using MD
simulations. It represents a mixture of the fluid and solute molecules as a LJ
fluid. The interaction parameters are determined by the law of the
corresponding states. Thus, knowing the critical parameters of the solvent and
solute, we can fit the critical point of the LJ fluid, obtaining the
parameters of the potential of the solute-solute and solvent-solvent
interactions. The parameters of the solute-solvent interaction potential are
determined by the standard Berthelot-Lorenz mixing rules. The reduced
quantities can be obtained in the standard way as follows:
$T^{\ast}=T/\varepsilon$ \- temperature, $\rho^{\ast}=\rho\sigma^{3}$ \-
density and $P^{\ast}=P\sigma^{3}/\varepsilon$ \- pressure. The reduced
critical parameters for the pure and untruncated LJ fluid were calculated by
Monte Carlo simulations in the grand canonical ensemble [62]:
$T_{c}^{*}=1.3120$, $\rho_{c}^{*}=0.316$ and $P_{c}^{*}=0.1279$. The obtained
parameters of the LJ potential for CO2 and three studied compounds are
presented in Table 2.

The simulation procedure is the same as the one described in the previous
section for AAMD, except for the change in the set of the scaling $\lambda$
parameters. As we consider only the LJ contribution for such simulations, the
parameters take the following values
$\\{\lambda_{LJ},\lambda_{C}\\}$:$\\{0.0,0.0\\}$, $\\{0.1,0.0\\}$,
$\\{0.2,0.0\\}$, $\\{0.3,0.0\\}$, $\\{0.4,0.0\\}$, $\\{0.5,0.0\\}$,
$\\{0.6,0.0\\}$, $\\{0.7,0.0\\}$, $\\{0.8,0.0\\}$, $\\{0.9,0.0\\}$,
$\\{1.0,0.0\\}$.

Table 2: Parameters of the LJ potential of the CO2 and drug compounds used in the CGMD simulations and obtained from the fitting of the LJ fluid liquid-vapor critical point. compound |  $\sigma$, $[$\mathrm{\SIUnitSymbolAngstrom}$]$ |  $\varepsilon$, $[K]$
---|---|---
CO2 | 3.673 | 231.805
aspirin | 6.789 | 581.479
ibuprofen | 7.565 | 571.387
carbamazepine | 7.439 | 599.718

## Experimental solvation free energy

To compare the calculated values with the experimental results we propose to
extract the solvation free energy values from the experimentally measured
solubility data available in literature. As the studied compounds are
sparingly soluble in a supercritical fluid, we consider the limit of their
infinite dilution. One should also assume that a supercritical fluid does not
dissolve in the solid phase, the solute molar volume does not depend on
pressure and the fugacity coefficient of the solid phase can be equal to unity
due to the fact that the sublimation pressure of the solute is rather small.
Thus, following from the condition of the chemical equilibrium for the solute
between its solid and solution phases (the details of derivation and
discussion can be found elsewhere [1, 38, 63, 64, 65]), one can obtain the
solvation free energy in the form

$\Delta G_{solv}=RT\ln{\frac{P^{sat}}{y_{2}\rho_{b}RT}}+v^{s}(P-P^{sat}),$
(10)

where $\Delta G_{solv}$ is the solvation free energy of the compound in the
scCO2 medium, $R$ is the gas constant, $T$ is the temperature, $P^{sat}$ is
the sublimation pressure of the solute, $y_{2}$ is the solubility of the
compound expressed in molar fraction, $\rho_{b}$ is the bulk density of the
fluid, $v^{s}$ is the molar volume of the solute and $P$ is the total
pressure, imposed in the system. It can be concluded from expression (10) that
apart from the experimental data on solubility, one has to know the properties
of the pure compound, i.e. the sublimation pressure and the molar volume. The
former can be estimated by semi-empirical methods (see, for example, [66, 67]
and cites therein) or, which is more preferable, can be measured
experimentally [68, 69, 70], and the latter one is often estimated by the
group contribution approaches (see [71]). The sources of the experimental data
discussed in the paper for all the solutes are presented in Table 3.

In Fig. 1 we show a comparison of the solvation free energy values obtained
from the ibuprofen (IBU) solubility experiment [72], using three sets of
sublimation pressure values, which we will designate as set Ardjmand_setI[72],
Ardjmand_setII [14] and Ardjmand_setIII [68] (Table 4). The error of the final
solvation free energy measurements was calculated as the uncertainty of the
indirect measurements, the main contribution to which turned out to mainly
depend on the uncertainty of the sublimation pressure values. As it can be
seen, the data are rather scattered and several values reach the discrepancies
of up to 1.5 $kcal/mol$. The molar volume was estimated [73] to be
$182.1\times 10^{-6}$ $m^{3}/mol$. It should be noted that the values of set
III were measured experimentally by the transpiration method [74], the values
of set I were numerically approximated [75], and in the case of set II the
sublimation pressure was considered as an adjustable parameter, the data for
which were obtained by minimizing the discrepancy between the experimental
solubility data and those correlated by the equation of state. In the case
when strict experimental data are available in literature, which is the case
for all the compounds studied here, we prefer to rely on them. Thus, in the
following sections we will utilize only the experimentally measured
sublimation pressure values.

Table 3: Experimental solvation free energy data. label | compound |  solubility data |  |
---|---|---|---|---
source |  sublimation pressure |  |  |
source |  molar volume |  |  |
source |  |  |  |
Ardjmand_setI | IBU | [72] | [72] | [72]
Ardjmand_setII | IBU | [72] | [14] | [72]
Ardjmand_setIII | IBU | [72] | [68] | [72]
Charoenchaitrakool | IBU | [76] | [69] | [72]
Huang | ASP | [77] | [69] | [77]
Champeau | ASP | [78] | [69] | [77]
Yamini | CBZ | [79] | [70] | [80]
IR | CBZ | [52] | [70] | [80]
Table 4: Sublimation pressure values of the IBU compound extracted from literature. |  $P^{sat}$, Pa
---|---
T, K | Ardjmand_setI [72] | Ardjmand_setII [14] | Ardjmand_setIII [68]
308.15 | 0.0495 | 0.0033 | 0.0081
313.15 | 0.0897 | 0.0075 | 0.0167
318.15 | 0.1600 | 0.0165 | 0.0335
Figure 1: Comparison of the IBU solvation free energy values, extracted from
the experiment [72], using three sets of the sublimation pressure data from
Table 4.

## Choice of the critical parameters

One of the advantages of the two proposed coarse-graining approaches is a
small number of input parameters. In fact, one only needs to know the solute
and solvent critical parameters for the system of choice. Finding them for a
solvent is usually not difficult as accurate information for the most
frequently used ones can be obtained from NIST [42]. The bottleneck of the
approach is the fact that the critical parameters of a solute cannot be
determined experimentally. It is common to estimate them by the group
contribution methods extensively used for this purpose [81, 82, 83]. However,
the values of the critical parameters are highly dependent on the calculation
method. Thus, the final solvation free energy values can vary a lot, as we
show in Fig.2, where we compare the IBU solvation free energy values obtained
by the cDFT approach for three isotherms using three sets of the estimated
critical parameters [76]. The symbols on the plot corresponding to the
experimental data suggest that the Ambrose method in this case gives the most
adequate critical values. Nevertheless, it is obvious that the final values of
the solvation free energy are highly dependent on the choice of the critical
parameters and the change in the critical temperature seems to affect the
outcomes more crucially than the change in the critical pressure. It has been
recently shown [84] that a possible solution to this ambiguity problem could
be the method of the critical parameters estimation proposed by G.Kontogeorgis
et al. [85] that more accurately predicts the sublimation enthalpy and vapor
pressure than the methods developed by K.Joback [82], K.Klincewicz [81] and
L.Constantinou [83]. Moreover, the approach by Kontogeorgis is applicable in
the cases of the compounds with pharmaceutical-like high molecular weights.
However, while the latter three methods use the value of the boiling point,
which in practice cannot be determined experimentally and thus has to be
estimated as well, Kontogeorgis’ method requires at least one experimental
vapor pressure data point.

Figure 2: Comparison of the IBU solvation free energy values obtained by the
cDFT approach using three sets of estimated critical parameters [76];
experimental data from [72, 76].

## RESULTS AND DISCUSSION

## Coexistence curves

Before turning to the results that we have obtained within the cDFT- and CGMD-
based approaches, it is instructive to understand how much our models of
solvent deviate from the real one. To understand this, we computed coexistence
curves using the equation of state implemented in our cDFT-based approach and
compared it to the data of the CO2 phase behavior extracted from NIST and to
those obtained by the LJ fluid critical point fitting procedure [62, 86]. As
one can see in the $\mathrm{T}$-$\mathrm{\rho}$ coordinates in Fig.3(a) we
have a slight deviation of the density at the chosen temperatures, with both
models overestimating the vapor branch and underestimating the liquid one. In
the case of the $\mathrm{P}$-$\mathrm{\rho}$ coexistence curves in Fig.3(b),
it can be seen that the data obtained by the equation of state overestimate
the critical pressure by over 20 bar, and the LJ fluid model overestimates it
by almost 10 bar. Thus, as we are interested in the supercritical region and
the model fluids do not describe the critical pressure accurately, it must be
said that the results obtained by the cDFT approach for the pressures lower
than approximately 95 bar and by the CGMD technique for the pressures lower
than 80 bar may not be reliable. It must be said that the complex molecular
potentials utilized in the AAMD simulations, in our case these are the Zhang
and TraPPE models, also demonstrate a distinct divergence from the
experimental coexistence curves [47, 87].

(a)

(b)

Figure 3: $T-\rho$ (3(a)) and $P-\rho$ (3(b)) coexisting curves of CO2
obtained from the EOS of the cDFT-based approach (black squares), extracted
from NIST (empty squares) and taken from the data of the LJ fluid [86] (red
squares). As it is seen, the model fluids overestimate the critical pressure
of the real solvent.

## Comparison of radial distribution functions

Another interesting point is the comparison of the structural properties of
model fluids. Fig.6 shows ASP-CO2 radial distribution functions (RDF) computed
by the cDFT approach and CGMD simulations with the corresponding potential
parameters (see Tables 1,2). The difference in the location of the first peaks
is approximately 0.2 $\mathrm{\SIUnitSymbolAngstrom}$, which basically follows
from the differences in the parameters of the potential. We also computed the
ASP-CO2 center of mass RDF by the AAMD simulations, using the OPLS-AA and GAFF
force-fields and Zhang model of CO2. Although the location of most peaks and
their heights differ drastically, it is interesting to note that the location
of the second peak of the OPLS-AA force field-based RDF is quite similar to
that of the first peak of the RDF obtained within the CGMD approach.
Considering that the solute molecule during the all-atom calculation is not
represented as a sphere, the first peak supposedly corresponds to the case
where the CO2 molecules approach the ASP molecule from the sides where the
effective distance between the molecules is shorter than the hard sphere
effective radius. In such a manner we can conclude that the structural
description of our model fluids resembles the all-atom center of mass MD
simulations results. As it will be seen from the following results, the
accurate description of the thermodynamic properties, namely, the liquid-vapor
coexistence curve, is more vital for the correct solvation free energy
evaluation than the exact account of the local fluid structure.

Figure 4: Qualitative comparison of the ASP-CO2 RDFs calculated by the coarse-
graining approaches (cDFT - red line, CGMD - black line) with the center-of-
mass RDFs computed by the AAMD simulations (force fields: OPLS-AA - blue line
and GAFF - green line) for T=308.15 K, P=150 bar.

## Uncertainties of the AAMD simulations

We would like to underline that the proposed coarse-graining techniques
suggest not only the reduction in the time and resources needed to compute the
solvation free energy values, but also reduce the ambiguity of the system
adjustment prior to the following calculation, as in the case of the AAMD
simulation. As we have pointed out above, the results of the solvation free
energy calculation using the AAMD simulations heavily rely on the choice of
the force field for the solute and solvent molecules [8, 1] and the method of
the partial atomic charges calculation [9]. In this section, we would like to
discuss the drawbacks of the AAMD approach we noticed while conducting the
simulation for the ASP compound. Fig. 5(a) shows the differences in the
solvation free energy values of the ASP molecule in scCO2 for three isotherms
using two force fields: GAFF and OPLS-AA and the Zhang model of CO2 [47]. As
one can see, the discrepancies between the values we have obtained with the
different force fields, are quite substantial, growing up to 2
$\mathrm{kcal/mol}$ at most. Fig. 5(b) shows the aspirin solvation free energy
values computed using the OPLS-AA force field and two different solvent
models: TraPPE and Zhang’s ones. The qualitative behavior is the same as shown
in paper [1], i.e. the values obtained with the TraPPE model are higher than
those obtained with the Zhang one, although the discrepancies are much larger.
Beside the significant differences in the results obtained using different
force fields and models of the system, it is important to note that the
absolute values of the solvation free energy we have obtained are quite high
compared with the experimental results, indicating rather low possibility of
achieving valid results without a thorough preliminary parametrization
routine.

(a)

(b)

Figure 5: Comparison of the ASP solvation free energy values for three
isotherms obtained from AAMD simulations with OPLS-AA (open circles) and GAFF
(filled circles) force fields and the Zhang model for CO2 (5(a)) and TraPPE
(filled circles) and Zhang (open circles) CO2 models with the OPLS-AA force-
field (5(b)); experimental data (crossed circles) from [77].

## CGMD simulation results

Let us now turn to the results of the proposed CGMD simulation technique.
Fig.6 shows a comparison of the solvation free energy values obtained within
the proposed CGMD approach with the interaction parameters presented in Table
2, AAMD simulations and the data extracted from the solubility experiments for
IBU at 308.15, 313.15 and 318.15 K (6(a)), carbamazepine (CBZ) at 308.15,
318.15 and 328.15 K (6(b)), and ASP at 308.15, 318.15 and 328.15 K (6(c)). The
experimental data were taken from the papers (see also table 3): IBU [76, 72],
CBZ [79] and ASP [77]. Regarding the AAMD simulations for the first solute we
used the Zhang model to describe the solvent and the OPLS force field for IBU,
which was fitted to quantum chemical simulation in scCO2 [88], thus leading to
much better agreement with the experimental data than the two others, where
there was no preliminary parameterization; for the second solute we used the
generic GROMOS 54A7 force-field [89] with the partial atomic charges of CBZ
computed by the Merz-Kollman method [59, 60], using the Gaussian 09 software
[61] with the PBE functional and 6-311++g(2d,p) basis set, which led to better
agreement with the experimental data; for the case of ASP we show here the
combination of the OPLS solute force field with Zhang’s CO2 model as the one
providing results that are closest to the experimental data. As one can see,
the results are rather surprising in the sense that a crude CGMD-based
approach can show reasonably good agreement with the experiment with the
largest deviation being near 1 kcal/mol, even outperforming the AAMD
procedures.

(a) IBU

(b) CBZ

(c) ASP

Figure 6: Comparison of the solvation free energy values obtained within the
CGMD approach (open circles), AAMD (open squares) and data extracted from the
solubility experiment (open triangles) for IBU, experiment data[76, 72]
(6(a)); CBZ, experiment data[79] (6(b)); ASP, experiment data[77] (6(c)).

## cDFT results

Now let us compare the results of the described cDFT approach with the results
of the experimental measurements and CGMD simulations. Fig.7 shows the values
of the solvation free energy of ASP in scCO2 obtained from the cDFT
calculations, CGMD and from the data of the solubility measurement experiments
[78, 77]. The ASP molar volume was estimated by the group-contribution method
and equals $129.64\times 10^{-6}$ $\mathrm{m^{3}/mol}$ [71]. One can observe
quite good agreement of the cDFT results with the experimental data,
especially in the region of the higher pressure values, while the CGMD results
are slightly understated.

Figure 7: Comparison of the ASP solvation free energy values obtained within
the cDFT approach (solid lines), CGMD (open circles) and experimental results
(open triangles) [78, 77].

Fig.8 shows the values of the carbamazepine solvation free energy obtained
within the cDFT approach, CGMD and experimental studies [79, 52]. Overall it
can be concluded that the agreement is quantitatively decent for the isotherms
at 308.15, 318.15 and 328.15 K. The comparison for the isotherms at
333.15-373.15 K does not demonstrate such sufficient agreement with the
experiment, although the higher the temperature and the pressure, the better
the agreement. Since the infrared (IR) spectroscopy experiment was conducted
under the isochoric conditions [52], we cannot propose a straightforward
comparison between the two experimental sets. Nevertheless, the quantitative
comparison of the solvation free energy values of several closely located
points from both sets shows that the IR data exceed Yamini’s data set
threefold at most.

Figure 8: Comparison of the CBZ solvation free energy values obtained within
the cDFT approach (solid lines), CGMD (open circles) and experimental results
(open triangles) [79, 52].

Fig. 9 demonstrates a comparison of the IBU solvation free energy values
obtained within the cDFT approach, CGMD simulations and experimental data [72,
76]. Here the critical parameters of IBU were estimated by the Ambrose method.
The correspondence with the experimental data is not as pronounced, although
it should be noted that the agreement between the cDFT results and the
experimental data for all the compounds is quite satisfactory at the pressures
from approximately 140 bar. For the case of IBU, the majority of the
experimental data points lie in the region below this pressure value. On the
other hand, when carefully analyzed, the experimental data for the highest
isotherm are rather scattered.

Figure 9: Comparison of the IBU solvation free energy values obtained within
the cDFT approach (solid lines), CGMD (open circles) and experimental results
(open triangles) [72, 76].

## CONCLUSION

We have proposed two coarse-graining approaches for fast computation of
solvation free energies of poorly soluble drug compounds in a scCO2 medium.
The first one is based on the utilization of the classical density functional
theory and the second one is built upon the use of MD simulations of a
supercritical solution represented as a Lennard Jones fluid. The
rationalization of such approximation lies in the fact that the main
contribution to the interactions between the solute and solvent molecules for
the described systems of sparingly soluble compounds in a scCO2 medium is that
of the van der Waals type. It seems more crucial to model the fluid phase
diagram correctly in this region by the fitting of the critical liquid-vapor
point than to make a detailed description of the internal structure of the
molecules. The satisfying agreement reached between the obtained results and
the experimental ones can be a reason to believe that such approaches can be
applied as an acceptable alternative to the all-atom MD simulations.

##  Appendix A Calculation of the solvation free energy via cDFT approach

To calculate the solvation Gibbs free energy of the solute in the
supercritical fluid, let us start from the expression for the grand
thermodynamic potential of scCO2 in the external field with the potential
energy $V_{ext}(\mathbf{r})$, created by the fixed solute molecule that is
posed at the origin

$\Omega[\rho(\mathbf{r})]=F_{int}[\rho(\mathbf{r})]+\int
d\mathbf{r}\rho(\mathbf{r})(V_{ext}(\mathbf{r})-\mu),$ (11)

where $F_{int}[\rho(\mathbf{r})]$ is the intrinsic Helmholtz free energy of
the scCO2 and $\mu$ is the chemical potential of bulk fluid at certain
temperature and pressure. The intrinsic Helmholtz free energy, in its turn,
can be comprised of two contributions

$F_{int}[\rho(\mathbf{r})]=k_{B}T\int
d\mathbf{r}\rho(\mathbf{r})[\ln(\Lambda^{3}\rho(\mathbf{r}))-1]+F_{ex}[\rho(\mathbf{r})],$
(12)

where the first one is the Helmholtz free energy of the ideal gas and the
second one is the excess Helmholtz free energy of fluid; $\Lambda$ is the
thermal de Broglie wavelength. We assume the CO2 molecules as coarse-grained
hard spheres interacting thorough the effective LJ pairwise potential

$V_{ff}(r)=4\varepsilon_{ff}\left(\left(\frac{\sigma_{ff}}{r}\right)^{12}-\left(\frac{\sigma_{ff}}{r}\right)^{6}\right)\Theta(r_{c}-r),$
(13)

where the $\varepsilon_{ff}$ and $\sigma_{ff}$ are LJ parameters of the
interaction between two molecules of the fluid and $\Theta(x)$ is the
Heaviside step function. The potential is divided into two parts at the
minimum of the LJ potential $r_{m}=2^{1/6}\sigma_{ff}$ in accordance with the
WCA procedure [90], and cut at $r_{c}=5\sigma_{ff}$. Thus, the excess part of
the free energy can be written as a sum of two contributions as follows

$F_{ex}[\rho(\mathbf{r})]=F_{hs}[\rho(\mathbf{r})]+F_{att}[\rho(\mathbf{r})],$
(14)

where the $F_{hs}[\rho(\mathbf{r})]$ is the Helmholtz free energy of the hard
sphere system and $F_{att}[\rho(\mathbf{r})]$ is the contribution of the
attractive interaction between the molecules of the fluid. In its turn, the
hard sphere part of the excess free energy is approximated with the help of
Rosenfeld’s version of the Fundamental Measure Theory (FMT) [32] in a certain
way

$F_{hs}[\rho(\mathbf{r})]=k_{B}T\int d\mathbf{r}\Phi(\mathbf{r}),$ (15)

where the $\Phi(\mathbf{r})$ is the excess free energy density, which is given
by the following expression

$\Phi=-n_{0}\ln(1-n_{3})+\frac{n_{1}n_{2}-{n_{V1}}\cdot{n_{V2}}}{1-n_{3}}+\frac{n_{2}^{3}-3n_{2}{n_{V1}}\cdot{n_{V2}}}{24\pi(1-n_{3})^{2}},$
(16)

where $n_{i}\cdot n_{j}$ denotes scalar multiplication.

The excess Helmholtz free energy density $\Phi(\mathbf{r})$ is a function of
the weighted densities

$n_{\alpha}(\mathbf{r})=\int
d\mathbf{r}^{\prime}\rho(\mathbf{r}^{\prime})\omega^{(\alpha)}(\mathbf{r}-\mathbf{r}^{\prime}),$
(17)

where the $\omega^{(\alpha)}(\mathbf{r}-\mathbf{r}^{\prime})$ are the weighted
functions, which characterize the hard sphere geometry: the volume, the
surface area and the mean radius of the curvature. There are six such
functions for the case of the three-dimensional hard sphere system, but only
three of them are independent [32]

$\displaystyle\omega^{(3)}(\mathbf{r})$ $\displaystyle=\Theta(R-r)$
$\displaystyle\omega^{(2)}(\mathbf{r})$ $\displaystyle=\delta(R-r)$ (18)
$\displaystyle{\omega}^{(V2)}(\mathbf{r})$
$\displaystyle=\frac{\mathbf{r}}{r}\delta(R-r).$

Other three are dependent on the described above functions in the following
way: $\omega^{(1)}=\omega^{(2)}/4\pi R$, $\omega^{(0)}=\omega^{(2)}/4\pi
R^{2}$, $\omega^{(V1)}=\omega^{(V2)}/4\pi R$, where $\delta(r)$ is the Dirac
delta-function and $R=d_{BH}/2$ is the effective hard sphere radius. The
effective Barker-Henderson (BH) diameter is determined by the following Pade
approximation [41]

$d_{BH}=\sigma_{ff}\frac{1.068\varepsilon_{ff}/k_{B}T+0.3837}{\varepsilon_{ff}/k_{B}T+0.4293}.$
(19)

The weighted functions give rise to the corresponding weighted densities
$n_{\alpha}(\mathbf{r})$

$\displaystyle n_{3}(r)$
$\displaystyle=\frac{\pi}{r}\int\limits_{r-R}^{r+R}dr^{\prime}r^{\prime}\rho(r^{\prime})[R^{2}-(r-r^{\prime})^{2}]$
$\displaystyle n_{2}(r)$ $\displaystyle=\frac{2\pi
R}{r}\int\limits_{r-R}^{r+R}dr^{\prime}r^{\prime}\rho(r^{\prime})$ (20)
$\displaystyle n_{V2}(r)$
$\displaystyle=\frac{\pi}{r^{2}}\int\limits_{r-R}^{r+R}dr^{\prime}r^{\prime}\rho(r^{\prime})[r^{2}+R^{2}-r^{\prime
2}]$

The contribution of the interparticle attraction to the excess free energy is
determined by the mean-field approximation

$F_{att}[\rho(\mathbf{r})]=\frac{1}{2}\int d\mathbf{r}\rho(\mathbf{r})\int
d\mathbf{r}^{\prime}\rho(\mathbf{r}^{\prime})\phi_{WCA}(\mathbf{r}-\mathbf{r}^{\prime}),$
(21)

where the effective WCA pair potential of attractive interactions is

$\phi_{WCA}(r)=\left\\{\begin{array}[]{lr}-\varepsilon_{ff},&r<r_{m}\\\
4\varepsilon_{ff}\left[\left(\frac{\sigma_{ff}}{r}\right)^{12}-\left(\frac{\sigma_{ff}}{r}\right)^{6}\right],&r_{m}<r<r_{c}.\end{array}\right.$
(22)

We model the solute molecule as a center generating external LJ potential

$V_{ext}(\mathbf{r})=4\varepsilon_{sf}\left(\left(\frac{\sigma_{sf}}{r}\right)^{12}-\left(\frac{\sigma_{sf}}{r}\right)^{6}\right)$
(23)

with the effective parameters of interaction between the solute and fluid
molecules. We obtain these parameters with the help of the Berthelot-Lorenz
mixing rules: $\sigma_{sf}=(\sigma_{ss}+\sigma_{ff})/2$ and
$\varepsilon_{sf}=\sqrt{\varepsilon_{ss}\varepsilon_{ff}}$, where the
parameters of the interaction between two molecules of the active compound
($\sigma_{ss}$, $\varepsilon_{ss}$) and CO2 ($\sigma_{ff}$,
$\varepsilon_{ff}$) can be obtained from the fitting of the respective
critical parameters111Critical density and temperature - for CO2, critical
pressure and temperature - for the compound. of the liquid-gas transition,
which values are available in the literature (see below).

For the bulk phase of the fluid, i.e. when $V_{ext}(\mathbf{r})=0$ and
$\rho(\mathbf{r})=\rho_{b}=const$, the functional (12) can be reduced to the
following form

$f=\frac{F_{int}}{V}=f_{id}+f_{ex},$ (24)

where $V$ is the system volume,
$f_{id}=\rho_{b}k_{B}T(\ln(\Lambda^{3}\rho_{b})-1)$ is the ideal free energy
density, and the excess contribution to the Helmholtz free energy density of
the bulk fluid takes the following form

$f_{ex}(\rho_{b},T)=\rho_{b}k_{B}T\left(-\ln(1-\eta)+\frac{3\eta}{1-\eta}+\frac{3\eta^{2}}{2(1-\eta)^{2}}\right)+\frac{1}{2}B_{WCA}\rho_{b}^{2},$
(25)

where $\eta=\pi d_{BH}^{3}\rho_{b}/6$ is the hard sphere system packing
fraction and the following auxiliary function, corresponding to the attractive
contribution, is introduced

$B_{WCA}=-\frac{32\sqrt{2}}{9}\pi\varepsilon_{ff}\sigma_{ff}^{3}+\frac{16}{3}\pi\varepsilon_{ff}\sigma_{ff}^{3}\left[\left(\frac{\sigma_{ff}}{r_{c}}\right)^{3}-\frac{1}{3}\left(\frac{\sigma_{ff}}{r_{c}}\right)^{9}\right].$
(26)

The total pressure then takes the following form

$P=\rho_{b}\frac{\partial
f}{\partial\rho_{b}}-f=\rho_{b}k_{B}T\frac{1+\eta+\eta^{2}}{(1-\eta)^{3}}+\frac{1}{2}B_{WCA}\rho_{b}^{2}.$
(27)

One can obtain LJ parameters of interaction between the fluid molecules as a
result of solving following system of the equations, using known critical
parameters of CO2

$\frac{\partial P}{\partial\rho_{b}}\Bigr{|}_{\begin{subarray}{c}T=T_{c}\\\
\rho_{b}=\rho_{c}\end{subarray}}=0,~{}~{}\frac{\partial^{2}P}{\partial\rho_{b}^{2}}\Bigr{|}_{\begin{subarray}{c}T=T_{c}\\\
\rho_{b}=\rho_{c}\end{subarray}}=0.$ (28)

In the same way the LJ parameters for the solute can be found.

In order to obtain the density profile, it is necessary to minimize the grand
thermodynamic potential with respect to $\rho(\mathbf{r})$ and solve
numerically the Euler-Lagrange equation

$\frac{\delta\Omega[\rho(\mathbf{r})]}{\delta\rho(\mathbf{r})}=0,$ (29)

this leads to the following expression

$\rho(\mathbf{r})=\rho_{b}\exp\left[\frac{\mu_{ex}(\rho_{b},T)-c^{(1)}_{fmt}(\mathbf{r})+\int
d\mathbf{r}^{\prime}\rho(\mathbf{r}^{\prime})\phi_{WCA}(\mathbf{r}-\mathbf{r}^{\prime})-V_{ext}(\mathbf{r})}{k_{B}T}\right],$
(30)

where $\mu_{ex}(\rho_{b},T)$ is the excess chemical potential of the bulk
phase, $c^{(1)}_{fmt}(\mathbf{r})=\delta
F_{hs}[\rho(\mathbf{r})]/{\delta\rho(\mathbf{r})}$ is the one-particle direct
correlation function of the hard-sphere system within the FMT. From the
definition of the one-particle direct correlation function and from the
obtained expression for the excess free energy within the FMT, we can write
the following

$c^{(1)}_{fmt}(\mathbf{r})=k_{B}T\int\limits_{r-R}^{r+R}dr^{\prime}\left[\frac{\pi}{r}\frac{\partial\Phi}{\partial
n_{3}}r^{\prime}(R^{2}-(r-r^{\prime})^{2})+\frac{2\pi
R}{r}\frac{\partial\Phi}{\partial
n_{2}}r^{\prime}+\frac{\partial\Phi}{\partial n_{V2}}\frac{\pi}{r}(r^{\prime
2}-r^{2}+R^{2})\right],$ (31)

where the derivatives are determined as follows

$\displaystyle\frac{\partial\Phi}{\partial n_{3}}$
$\displaystyle=\frac{n_{2}}{4\pi
R^{2}(1-n_{3})}+\frac{n_{2}^{2}-n_{V2}^{2}}{4\pi
R(1-n_{3})^{2}}+\frac{n_{2}^{3}-3n_{2}n_{V2}^{2}}{12\pi(1-n_{3})^{3}},$
$\displaystyle\frac{\partial\Phi}{\partial n_{2}}$
$\displaystyle=-\ln(1-n_{3})/(4\pi R^{2})+\frac{n_{2}}{2\pi
R(1-n_{3})}+\frac{n_{2}^{2}-n_{V2}^{2}}{8\pi(1-n_{3})^{2}},$ (32)
$\displaystyle\frac{\partial\Phi}{\partial n_{V2}}$
$\displaystyle=-\frac{n_{V2}}{2\pi
R(1-n_{3})}-\frac{n_{2}n_{V2}}{4\pi(1-n_{3})^{2}}.$

In order to calculate the contribution of the attractive interactions between
two fluid particles in spherical geometry, one have to distinguish the
situations where $r-r_{c}<0$ and $r-r_{c}>0$. Also, using the spherical
symmetry of the system, we can calculate the integral only for
$\mathbf{r}=r\mathbf{z}$ without loss of generality. Here we will present the
equations for the latter case as an example

$\displaystyle\int
d\mathbf{r}^{\prime}\rho(\mathbf{r}^{\prime})\phi_{WCA}(\mathbf{r}-\mathbf{r}^{\prime})$
$\displaystyle=$ $\displaystyle
2\pi\varepsilon_{ff}\int_{r-r_{c}}^{r-r_{m}}d{r}^{\prime}{r}^{\prime}\rho({r}^{\prime})\frac{1}{r}\biggl{[}\frac{2\sigma_{ff}^{12}}{5(r-{r}^{\prime})^{10}}-\frac{\sigma_{ff}^{6}}{(r-{r}^{\prime})^{4}}-\frac{2\sigma_{ff}^{12}}{5r_{c}^{10}}+\frac{\sigma_{ff}^{6}}{r_{c}^{4}}\biggr{]}$
(33) $\displaystyle+$ $\displaystyle
2\pi\varepsilon_{ff}\int_{r+r_{m}}^{r+r_{c}}d{r}^{\prime}{r}^{\prime}\rho({r}^{\prime})\frac{1}{r}\biggl{[}\frac{2\sigma_{ff}^{12}}{5(r-{r}^{\prime})^{10}}-\frac{\sigma_{ff}^{6}}{(r-{r}^{\prime})^{4}}-\frac{2\sigma_{ff}^{12}}{5r_{c}^{10}}+\frac{\sigma_{ff}^{6}}{r_{c}^{4}}\biggr{]}$
$\displaystyle-$
$\displaystyle\frac{\pi\varepsilon_{ff}}{r}\int_{r-r_{m}}^{r+r_{m}}d{r}^{\prime}{r}^{\prime}\rho({r}^{\prime})[r_{m}^{2}-(r-{r}^{\prime})^{2}]$
$\displaystyle+$ $\displaystyle
2\pi\varepsilon_{ff}\int_{r-r_{m}}^{r+r_{m}}d{r}^{\prime}{r}^{\prime}\rho({r}^{\prime})\frac{1}{r}\biggl{[}\frac{2\sigma_{ff}^{12}}{5r_{m}^{10}}-\frac{\sigma_{ff}^{6}}{r_{m}^{4}}-\frac{2\sigma_{ff}^{12}}{5r_{c}^{10}}+\frac{\sigma_{ff}^{6}}{r_{c}^{4}}\biggr{]}.$

Finally, the solvation Gibbs free energy can be calculated as the excess grand
thermodynamic potential as follows

$\Delta G_{solv}=\Omega[\rho(\mathbf{r})]-\Omega[\rho_{b}].$ (34)

##  Appendix B Aspirin partial atomic charges

Two most stable ASP conformers are represented at the Fig. 10. The partial
atomic charges, averaged over these two conformers, are shown in Table 5.

(a)

(b)

Figure 10: Two most stable ASP conformers as calculated from quantum chemistry. Table 5: Atomic partial charges of the ASP conformers, computed with the Merz-Kollman method, PBE functional, 6-311++g(2d,p) basis set. conf.1 | conf.2
---|---
atom | charge | atom | charge
H1 | 0.402490 | H1 | 0.394904
O2 | -0.551911 | O2 | -0.530321
C3 | 0.632537 | C3 | 0.613256
O4 | -0.513812 | O4 | -0.503629
C5 | -0.106338 | C5 | -0.118470
C6 | 0.374670 | C6 | 0.396578
O7 | -0.416109 | O7 | -0.418618
C8 | 0.768870 | C8 | 0.772999
O9 | -0.502975 | O9 | -0.500559
C10 | -0.557154 | C10 | -0.600964
H11 | 0.152293 | H11 | 0.166754
H12 | 0.194728 | H12 | 0.192129
H13 | 0.155861 | H13 | 0.172574
C14 | -0.259618 | C14 | -0.278193
H15 | 0.165669 | H15 | 0.169621
C16 | -0.081902 | C16 | -0.083413
H17 | 0.136997 | H17 | 0.139358
C18 | -0.156598 | C18 | -0.161393
H19 | 0.136687 | H19 | 0.139205
C20 | -0.137291 | C20 | -0.116117
H21 | 0.162907 | H21 | 0.154300

##  Appendix C Experimental solvation free energy data

See pages - of table_for_si.pdf

### ACKNOWLEDGMENTS

The research was funded by the Ministry of Science and Higher Education of the
Russian Federation (grant no. RFMEFI61618$\times$0097). This research was done
on the supercomputer facilities provided by NRU HSE.

## References

  * [1] J. Noroozi, C. Ghotbi, J. J. Sardroodi, J. Karimi-Sabet, and M. A. Robert, “Solvation free energy and solubility of acetaminophen and ibuprofen in supercritical carbon dioxide: Impact of the solvent model,” The Journal of Supercritical Fluids, vol. 109, pp. 166–176, 2016.
  * [2] A. I. Frolov, “Accurate calculation of solvation free energies in supercritical fluids by fully atomistic simulations: Probing the theory of solutions in energy representation,” Journal of chemical theory and computation, vol. 11, no. 5, pp. 2245–2256, 2015.
  * [3] S. Bruckner and S. Boresch, “Efficiency of alchemical free energy simulations. ii. improvements for thermodynamic integration,” Journal of computational chemistry, vol. 32, no. 7, pp. 1320–1333, 2011.
  * [4] N. Hansen and W. F. Van Gunsteren, “Practical aspects of free-energy calculations: a review,” Journal of chemical theory and computation, vol. 10, no. 7, pp. 2632–2647, 2014.
  * [5] X. Jia, M. Wang, Y. Shao, G. König, B. R. Brooks, J. Z. Zhang, and Y. Mei, “Calculations of solvation free energy through energy reweighting from molecular mechanics to quantum mechanics,” Journal of chemical theory and computation, vol. 12, no. 2, pp. 499–511, 2016.
  * [6] M. Misin, M. V. Fedorov, and D. S. Palmer, “Hydration free energies of molecular ions from theory and simulation,” The Journal of Physical Chemistry B, vol. 120, no. 5, pp. 975–983, 2016.
  * [7] M. R. Shirts, “Best practices in free energy calculations for drug design,” in Computational drug discovery and design, pp. 425–467, Springer, 2012\.
  * [8] M. Lundborg and E. Lindahl, “Automatic gromacs topology generation and comparisons of force fields for solvation free energy calculations,” The Journal of Physical Chemistry B, vol. 119, no. 3, pp. 810–823, 2015.
  * [9] J. P. Jämbeck, F. Mocci, A. P. Lyubartsev, and A. Laaksonen, “Partial atomic charges and their impact on the free energy of solvation,” Journal of computational chemistry, vol. 34, no. 3, pp. 187–197, 2013.
  * [10] G. C. da Silva, G. M. Silva, F. W. Tavares, F. P. Fleming, and B. A. Horta, “Are all-atom any better than united-atom force fields for the description of liquid properties of alkanes?,” Journal of Molecular Modeling, vol. 26, no. 11, pp. 1–17, 2020.
  * [11] A. D. Glova, I. V. Volgin, V. M. Nazarychev, S. V. Larin, S. V. Lyulin, and A. A. Gurtovenko, “Toward realistic computer modeling of paraffin-based composite materials: critical assessment of atomic-scale models of paraffins,” RSC Advances, vol. 9, no. 66, pp. 38834–38847, 2019.
  * [12] K. D. Papavasileiou, L. D. Peristeras, A. Bick, and I. G. Economou, “Molecular dynamics simulation of pure n-alkanes and their mixtures at elevated temperatures using atomistic and coarse-grained force fields,” The Journal of Physical Chemistry B, vol. 123, no. 29, pp. 6229–6243, 2019.
  * [13] J. P. Ewen, C. Gattinoni, F. M. Thakkar, N. Morgan, H. A. Spikes, and D. Dini, “A comparison of classical force-fields for molecular dynamics simulations of lubricants,” Materials, vol. 9, no. 8, p. 651, 2016.
  * [14] C. Garlapati and G. Madras, “Temperature independent mixing rules to correlate the solubilities of antibiotics and anti-inflammatory drugs in scco2,” Thermochimica Acta, vol. 496, no. 1-2, pp. 54–58, 2009.
  * [15] E. Moine, R. Privat, J.-N. Jaubert, B. Sirjean, N. Novak, E. Voutsas, and C. Boukouvalas, “Can we safely predict solvation gibbs energies of pure and mixed solutes with a cubic equation of state?,” Pure and Applied Chemistry, vol. 91, no. 8, pp. 1295–1307, 2019.
  * [16] G. M. Kontogeorgis, X. Liang, A. Arya, and I. Tsivintzelis, “Equations of state in three centuries. are we closer to arriving to a single model for all applications?,” Chemical Engineering Science: X, vol. 7, p. 100060, 2020\.
  * [17] P. Hutacharoen, S. Dufal, V. Papaioannou, R. M. Shanker, C. S. Adjiman, G. Jackson, and A. Galindo, “Predicting the solvation of organic compounds in aqueous environments: from alkanes and alcohols to pharmaceuticals,” Industrial & Engineering Chemistry Research, vol. 56, no. 38, pp. 10856–10876, 2017.
  * [18] H. A. A. El, C. Si-Moussa, S. Hanini, and M. Laidi, “Application of pc-saft and cubic equations of state for the correlation of solubility of some pharmaceutical and statin drugs in sc-co2,” Chemical Industry and Chemical Engineering Quarterly/CICEQ, vol. 19, no. 3, pp. 449–460, 2013.
  * [19] M. H. Anvari and G. Pazuki, “A study on the predictive capability of the saft-vr equation of state for solubility of solids in supercritical co2,” The Journal of Supercritical Fluids, vol. 90, pp. 73–83, 2014.
  * [20] H. Yang and C. Zhong, “Modeling of the solubility of aromatic compounds in supercritical carbon dioxide–cosolvent systems using saft equation of state,” The Journal of supercritical fluids, vol. 33, no. 2, pp. 99–106, 2005.
  * [21] G. Sodeifian, S. A. Sajadian, and R. Derakhsheshpour, “Experimental measurement and thermodynamic modeling of lansoprazole solubility in supercritical carbon dioxide: Application of saft-vr eos,” Fluid Phase Equilibria, vol. 507, p. 112422, 2020.
  * [22] S. Z. Mahmoudabadi and G. Pazuki, “Application of pc-saft eos for pharmaceuticals: Solubility, co-crystal, and thermodynamic modeling,” Journal of Pharmaceutical Sciences, 2021.
  * [23] N. Ramírez-Vélez, A. Piña Martinez, J.-N. Jaubert, and R. Privat, “Parameterization of saft models: Analysis of different parameter estimation strategies and application to the development of a comprehensive database of pc-saft molecular parameters,” Journal of Chemical & Engineering Data, vol. 65, no. 12, pp. 5920–5932, 2020.
  * [24] Y. Shimoyama and Y. Iwai, “Development of activity coefficient model based on cosmo method for prediction of solubilities of solid solutes in supercritical carbon dioxide,” The Journal of Supercritical Fluids, vol. 50, no. 3, pp. 210–217, 2009.
  * [25] L.-H. Wang and S.-T. Lin, “A predictive method for the solubility of drug in supercritical carbon dioxide,” The Journal of Supercritical Fluids, vol. 85, pp. 81–88, 2014.
  * [26] A. V. Marenich, C. J. Cramer, and D. G. Truhlar, “Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions,” The Journal of Physical Chemistry B, vol. 113, no. 18, pp. 6378–6396, 2009.
  * [27] A. C. Chamberlin, D. G. Levitt, C. J. Cramer, and D. G. Truhlar, “Modeling free energies of solvation in olive oil,” Molecular pharmaceutics, vol. 5, no. 6, pp. 1064–1079, 2008.
  * [28] A. Klamt, COSMO-RS: from quantum chemistry to fluid phase thermodynamics and drug design.  Elsevier, 2005.
  * [29] M. Misin, D. S. Palmer, and M. V. Fedorov, “Predicting solvation free energies using parameter-free solvent models,” The Journal of Physical Chemistry B, vol. 120, no. 25, pp. 5724–5731, 2016.
  * [30] V. F. Sokolov and G. N. Chuev, “Fundamental measure theory of hydrated hydrocarbons,” Journal of molecular modeling, vol. 13, no. 2, pp. 319–326, 2007.
  * [31] G. Chuev and V. Sokolov, “Hydration of hydrophobic solutes treated by the fundamental measure approach,” The Journal of Physical Chemistry B, vol. 110, no. 37, pp. 18496–18503, 2006.
  * [32] Y. Rosenfeld, “Free-energy model for the inhomogeneous hard-sphere fluid mixture and density-functional theory of freezing,” Physical review letters, vol. 63, no. 9, p. 980, 1989.
  * [33] D.-Y. Peng and D. B. Robinson, “A new two-constant equation of state,” Industrial & Engineering Chemistry Fundamentals, vol. 15, no. 1, pp. 59–64, 1976.
  * [34] S. Zhao, R. Ramirez, R. Vuilleumier, and D. Borgis, “Molecular density functional theory of solvation: From polar solvents to water,” The Journal of chemical physics, vol. 134, no. 19, p. 194102, 2011.
  * [35] S. Zhao, Z. Jin, and J. Wu, “New theoretical method for rapid prediction of solvation free energy in water,” The Journal of Physical Chemistry B, vol. 115, no. 21, pp. 6971–6975, 2011.
  * [36] V. P. Sergiievskyi, G. Jeanmairet, M. Levesque, and D. Borgis, “Fast computation of solvation free energies with molecular density functional theory: Thermodynamic-ensemble partial molar volume corrections,” The journal of physical chemistry letters, vol. 5, no. 11, pp. 1935–1942, 2014.
  * [37] L. Gendre, R. Ramirez, and D. Borgis, “Classical density functional theory of solvation in molecular solvents: Angular grid implementation,” Chemical Physics Letters, vol. 474, no. 4-6, pp. 366–370, 2009.
  * [38] M. Baghbanbashi, N. Hadidi, and G. Pazuki, “Solubility of pharmaceutical compounds in supercritical carbon dioxide: Application, experimental, and mathematical modeling,” in Green Sustainable Process for Chemical and Environmental Engineering and Science, pp. 185–254, Elsevier, 2020.
  * [39] L. Padrela, M. A. Rodrigues, A. Duarte, A. M. Dias, M. E. Braga, and H. C. de Sousa, “Supercritical carbon dioxide-based technologies for the production of drug nanoparticles/nanocrystals–a comprehensive review,” Advanced drug delivery reviews, vol. 131, pp. 22–78, 2018.
  * [40] A. J. Archer, B. Chacko, and R. Evans, “The standard mean-field treatment of inter-particle attraction in classical dft is better than one might expect,” The Journal of chemical physics, vol. 147, no. 3, p. 034501, 2017.
  * [41] L. Verlet and J.-J. Weis, “Equilibrium theory of simple liquids,” Physical Review A, vol. 5, no. 2, p. 939, 1972.
  * [42] E. W. Lemmon, M. O. McLinden, and D. G. Friend, Thermophysical properties of fluid systems.  NIST chemistry WebBook, 1998.
  * [43] A. K. Malde, L. Zuo, M. Breeze, M. Stroet, D. Poger, P. C. Nair, C. Oostenbrink, and A. E. Mark, “An automated force field topology builder (atb) and repository: version 1.0,” Journal of chemical theory and computation, vol. 7, no. 12, pp. 4026–4037, 2011.
  * [44] M. Stroet, B. Caron, K. M. Visscher, D. P. Geerke, A. K. Malde, and A. E. Mark, “Automated topology builder version 3.0: prediction of solvation free enthalpies in water and hexane,” Journal of chemical theory and computation, vol. 14, no. 11, pp. 5834–5845, 2018.
  * [45] J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, and D. A. Case, “Development and testing of a general amber force field,” Journal of computational chemistry, vol. 25, no. 9, pp. 1157–1174, 2004.
  * [46] G. A. Kaminski, R. A. Friesner, J. Tirado-Rives, and W. L. Jorgensen, “Evaluation and reparametrization of the opls-aa force field for proteins via comparison with accurate quantum chemical calculations on peptides,” The Journal of Physical Chemistry B, vol. 105, no. 28, pp. 6474–6487, 2001\.
  * [47] Z. Zhang and Z. Duan, “An optimized molecular potential for carbon dioxide,” The Journal of chemical physics, vol. 122, no. 21, p. 214507, 2005.
  * [48] J. J. Potoff and J. I. Siepmann, “Vapor–liquid equilibria of mixtures containing alkanes, carbon dioxide, and nitrogen,” AIChE journal, vol. 47, no. 7, pp. 1676–1682, 2001.
  * [49] S. Canzar, M. El-Kebir, R. Pool, K. Elbassioni, A. K. Malde, A. E. Mark, D. P. Geerke, L. Stougie, and G. W. Klau, “Charge group partitioning in biomolecular simulation,” Journal of Computational Biology, vol. 20, no. 3, pp. 188–198, 2013.
  * [50] K. B. Koziara, M. Stroet, A. K. Malde, and A. E. Mark, “Testing and validation of the automated topology builder (atb) version 2.0: prediction of hydration free enthalpies,” Journal of computer-aided molecular design, vol. 28, no. 3, pp. 221–233, 2014.
  * [51] Y. Budkov, A. Kolesnikov, D. Ivlev, N. Kalikin, and M. Kiselev, “Possibility of pressure crossover prediction by classical dft for sparingly dissolved compounds in scco2,” Journal of Molecular Liquids, vol. 276, pp. 801–805, 2019.
  * [52] N. Kalikin, M. Kurskaya, D. Ivlev, M. Krestyaninov, R. Oparin, A. Kolesnikov, Y. Budkov, A. Idrissi, and M. Kiselev, “Carbamazepine solubility in supercritical co2: A comprehensive study,” Journal of Molecular Liquids, p. 113104, 2020.
  * [53] S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, et al., “Gromacs 4.5: a high-throughput and highly parallel open source molecular simulation toolkit,” Bioinformatics, vol. 29, no. 7, pp. 845–854, 2013.
  * [54] M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, “Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,” SoftwareX, vol. 1, pp. 19–25, 2015.
  * [55] H. Bekker, H. Berendsen, E. Dijkstra, S. Achterop, R. Vondrumen, D. Vanderspoel, A. Sijbers, H. Keegstra, and M. Renardus, “Gromacs-a parallel computer for molecular-dynamics simulations,” in 4th International Conference on Computational Physics (PC 92), pp. 252–256, World Scientific Publishing, 1993.
  * [56] H. J. Berendsen, D. van der Spoel, and R. van Drunen, “Gromacs: a message-passing parallel molecular dynamics implementation,” Computer physics communications, vol. 91, no. 1-3, pp. 43–56, 1995.
  * [57] L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, “Packmol: a package for building initial configurations for molecular dynamics simulations,” Journal of computational chemistry, vol. 30, no. 13, pp. 2157–2164, 2009.
  * [58] C. H. Bennett, “Efficient estimation of free energy differences from monte carlo data,” Journal of Computational Physics, vol. 22, no. 2, pp. 245–268, 1976.
  * [59] U. C. Singh and P. A. Kollman, “An approach to computing electrostatic charges for molecules,” Journal of Computational Chemistry, vol. 5, no. 2, pp. 129–145, 1984.
  * [60] B. H. Besler, K. M. Merz Jr, and P. A. Kollman, “Atomic charges derived from semiempirical methods,” Journal of Computational Chemistry, vol. 11, no. 4, pp. 431–439, 1990.
  * [61] M. Frisch, G. Trucks, H. Schlegel, G. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. Petersson, et al., “Gaussian 09, revision b. 01,” 2013.
  * [62] J. J. Potoff and A. Z. Panagiotopoulos, “Critical point and phase behavior of the pure fluid and a lennard-jones mixture,” The Journal of chemical physics, vol. 109, no. 24, pp. 10914–10920, 1998.
  * [63] R. Hartono, G. A. Mansoori, and A. Suwono, “Prediction of solubility of biomolecules in supercritical solvents,” Chemical Engineering Science, vol. 56, no. 24, pp. 6949–6958, 2001.
  * [64] S. V. de Melo, G. M. N. Costa, A. Viana, and F. Pessoa, “Solid pure component property effects on modeling upper crossover pressure for supercritical fluid process synthesis: A case study for the separation of annatto pigments using sc-co2,” The Journal of Supercritical Fluids, vol. 49, no. 1, pp. 1–8, 2009.
  * [65] Z. Su and M. Maroncelli, “Simulations of solvation free energies and solubilities in supercritical solvents,” The Journal of chemical physics, vol. 124, no. 16, p. 164506, 2006.
  * [66] A. Komkoua Mbienda, C. Tchawoua, D. Vondou, and F. Mkankam Kamga, “Evaluation of vapor pressure estimation methods for use in simulating the dynamic of atmospheric organic aerosols,” International Journal of Geophysics, vol. 2013, 2013.
  * [67] S. O’Meara, A. M. Booth, M. H. Barley, D. Topping, and G. McFiggans, “An assessment of vapour pressure estimation methods,” Physical Chemistry Chemical Physics, vol. 16, no. 36, pp. 19453–19469, 2014.
  * [68] G. L. Perlovich, S. V. Kurkov, L. K. Hansen, and A. Bauer-Brandl, “Thermodynamics of sublimation, crystal lattice energies, and crystal structures of racemates and enantiomers:(+)-and ($\pm$)-ibuprofen,” Journal of pharmaceutical sciences, vol. 93, no. 3, pp. 654–666, 2004.
  * [69] G. L. Perlovich, S. V. Kurkov, A. N. Kinchin, and A. Bauer-Brandl, “Solvation and hydration characteristics of ibuprofen and acetylsalicylic acid,” Aaps Pharmsci, vol. 6, no. 1, pp. 22–30, 2004.
  * [70] K. V. Drozd, A. N. Manin, A. V. Churakov, and G. L. Perlovich, “Novel drug–drug cocrystals of carbamazepine with para-aminosalicylic acid: Screening, crystal structures and comparative study of carbamazepine cocrystal formation thermodynamics,” CrystEngComm, vol. 19, no. 30, pp. 4273–4286, 2017.
  * [71] X. Cao, N. Leyva, S. R. Anderson, and B. C. Hancock, “Use of prediction methods to estimate true density of active pharmaceutical ingredients,” International journal of pharmaceutics, vol. 355, no. 1-2, pp. 231–237, 2008\.
  * [72] M. Ardjmand, M. Mirzajanzadeh, and F. Zabihi, “Measurement and correlation of solid drugs solubility in supercritical systems,” Chinese Journal of Chemical Engineering, vol. 22, no. 5, pp. 549–558, 2014.
  * [73] E. Baum, Chemical property estimation: theory and application.  CRC Press, 1997.
  * [74] X. Zielenkiewicz, G. Perlovich, and M. Wszelaka-Rylik, “The vapour pressure and the enthalpy of sublimation: determination by inert gas flow method,” Journal of thermal analysis and calorimetry, vol. 57, no. 1, pp. 225–234, 1999.
  * [75] W. J. Lyman, W. F. Reehl, and D. H. Rosenblatt, “Handbook of chemical property estimation methods,” 1990.
  * [76] M. Charoenchaitrakool, F. Dehghani, N. Foster, and H. Chan, “Micronization by rapid expansion of supercritical solutions to enhance the dissolution rates of poorly water-soluble pharmaceuticals,” Industrial & engineering chemistry research, vol. 39, no. 12, pp. 4794–4802, 2000.
  * [77] Z. Huang, W. D. Lu, S. Kawi, and Y. C. Chiew, “Solubility of aspirin in supercritical carbon dioxide with and without acetone,” Journal of Chemical & Engineering Data, vol. 49, no. 5, pp. 1323–1327, 2004.
  * [78] M. Champeau, J.-M. Thomassin, C. Jérôme, and T. Tassaing, “Solubility and speciation of ketoprofen and aspirin in supercritical co2 by infrared spectroscopy,” Journal of Chemical & Engineering Data, vol. 61, no. 2, pp. 968–978, 2016.
  * [79] Y. Yamini, J. Hassan, and S. Haghgo, “Solubilities of some nitrogen-containing drugs in supercritical carbon dioxide,” Journal of Chemical & Engineering Data, vol. 46, no. 2, pp. 451–455, 2001.
  * [80] J.-h. Li, Z. Huang, J.-l. Wei, and L. Xu, “A new optimization method for parameter determination in modeling solid solubility in supercritical co2,” Fluid Phase Equilibria, vol. 344, pp. 117–124, 2013.
  * [81] K. Klincewicz and R. Reid, “Estimation of critical properties with group contribution methods,” AIChE Journal, vol. 30, no. 1, pp. 137–142, 1984\.
  * [82] K. G. Joback and R. C. Reid, “Estimation of pure-component properties from group-contributions,” Chemical Engineering Communications, vol. 57, no. 1-6, pp. 233–243, 1987.
  * [83] L. Constantinou and R. Gani, “New group contribution method for estimating properties of pure compounds,” AIChE Journal, vol. 40, no. 10, pp. 1697–1710, 1994.
  * [84] S. A. Jahromi and A. Roosta, “Estimation of critical point, vapor pressure and heat of sublimation of pharmaceuticals and their solubility in supercritical carbon dioxide,” Fluid Phase Equilibria, vol. 488, pp. 1–8, 2019.
  * [85] G. M. Kontogeorgis, I. Smirlis, I. V. Yakoumis, V. Harismiadis, and D. P. Tassios, “Method for estimating critical properties of heavy compounds suitable for cubic equations of state and its application to the prediction of vapor pressures,” Industrial & engineering chemistry research, vol. 36, no. 9, pp. 4008–4012, 1997.
  * [86] A. Z. Panagiotopoulos, N. Quirke, M. Stapleton, and D. Tildesley, “Phase equilibria by simulation in the gibbs ensemble: alternative derivation, generalization and application to mixture and membrane equilibria,” Molecular Physics, vol. 63, no. 4, pp. 527–545, 1988.
  * [87] T. Merker, J. Vrabec, and H. Hasse, “Comment on “an optimized potential for carbon dioxide”[j. chem. phys. 122, 214507 (2005)],” The Journal of chemical physics, vol. 129, no. 8, p. 214507, 2008.
  * [88] I. Fedorova, D. Ivlev, and M. Kiselev, “Conformational lability of ibuprofen in supercritical carbon dioxide,” Russian Journal of Physical Chemistry B, vol. 10, no. 7, pp. 1153–1162, 2016.
  * [89] N. Schmid, A. P. Eichenberger, A. Choutko, S. Riniker, M. Winger, A. E. Mark, and W. F. van Gunsteren, “Definition and testing of the gromos force-field versions 54a7 and 54b7,” European biophysics journal, vol. 40, no. 7, pp. 843–856, 2011.
  * [90] H. C. Andersen, J. D. Weeks, and D. Chandler, “Relationship between the hard-sphere fluid and fluids with realistic repulsive forces,” Physical Review A, vol. 4, no. 4, p. 1597, 1971.