The equilibrium between liquid phases, or between liquid and vapor, is critical to understanding a wide range of chemical and biological phenomena; indeed, understanding and predicting such equilibria is essential for optimizing chemical products and processes [1]. Researchers have generally employed these types of equilibria to theoretically estimate a range of properties, including the aqueous solvation energy of polar and nonpolar species, the vapor pressure and boiling point of pure liquids, the equilibrium constant of carboxylic acids in the aqueous phase, and the partition coefficient of a solute between two immiscible liquids [2-8].

The partition coefficient of a solute between two immiscible liquids is one of the most significant properties used in this context. This property is particularly valuable for describing the behavior of amphiphilic molecules and species used in the pharmaceutical industry. The distribution of a molecule between water and a lipid material is typically described using the 1-octanol/water partition coefficient Ko/w, commonly expressed as ( $K_{o/w}=\log P$ ), which is regarded as the most important physicochemical property for describing the affinity of chemical species for a given solvent [9-11].

1-Octanol is widely accepted as the material whose behavior most closely resembles that of lipids found in nature, and the partition coefficient is the most widely used descriptor in structure-property relationship (QSPR) studies [12]. Concurrently, QSPR models emerge from a process that begins with the description of molecular structures and culminates in the prediction of molecular behavior across different phases [13-16]. Consequently, QSPR models have been developed using properties calculated by quantum mechanical methods and molecular dynamics [17-20].

A variety of quantum mechanical solvation models, including C-PCM (Conductor-Like Polarizable Continuum Model) and COSMO (Conductor-Like Screening Model), have been used to investigate and calculate the properties of chemical compounds in solvents [21-27]. In these methods, the solvent is represented by a continuous dielectric model, which provides a satisfactory description of molecular properties at low computational cost and yields favorable results. Numerous methods have also been developed to predict the molecular properties of solutes in different solvent media, including the molecular volume of liquids, vapor pressure, and solvation energies [28-32]. Accordingly, this study employs the C-PCM solvation method to calculate the structural and energetic properties of a set of organic molecules. Two quantitative structure-property relationship models were then obtained to estimate the 1-octanol/water partition coefficient ( $K_{o/w}=\log P)$ and the free energies of solvation for these groups of organic molecules. While solvation free energies are essential for understanding solute-solvent interactions, the partition coefficient is inherently a transfer property between two phases. The distinction between transfer and solvation free energies is critical because the former accounts for the competitive solvation of the solute by two immiscible solvents, including effects of activity coefficients, conformational changes, and specific solvent-solvent interactions that are not captured by single-solvent continuum models.

The solvation free energy ( $\mathsf{\Delta }{G}_{solv}$ ) is the free energy change required to transfer a molecule from the vacuum into the solvent [23] and can be defined by Equation (1):

The term $\mathsf{\Delta }{G}_{elec}$ is the electrostatic component arising from the interaction between the solute and the solvent medium. The variables $\mathsf{\Delta }{G}_{rep}$ and $\mathsf{\Delta }{G}_{disp}$ are the free energies of repulsion and dispersion between the solute and solvent. Finally, the term $\mathsf{\Delta }{G}_{cav}$ is the free energy required to construct the cavity of the solute in the solvent; the last two terms are negligible.

2.1. Electrostatic Free Energy

Under the influence of a solvent described by a dielectric continuum model, the Hamiltonian of the solute can be defined by Equation (2) as:

Here, ${\widehat{H}}^o$ is the Hamiltonian of the isolated solute, and the operator $\widehat{V}$ describes the solute-solvent interaction, which depends linearly on the state function $\Psi$ of the solute. In this case, the electrostatic free energy $\mathsf{\Delta }{G}_{elec}$ is determined by variational minimization of the free energy of the solute, given by Equation (3):

In general, the interaction operator $\widehat{V}$ is written in terms of apparent polarization charges. In each tessera $i$ , a charge $q_i$ appears, according to the conductor-like boundary condition:

where $V\left(\overrightarrow{r}\right)$ and $V_{q_i}\left(\overrightarrow{r}\right)$ are the electrostatic potentials due to the solute and to the polarization charges, respectively, and $\overrightarrow{r}$ is a point on the segment surface.

The vector of polarization charges $Q$ can be defined by Equation (5):

Here, the vector $V$ contains the electrostatic potential due to the solute on the surface segment. The elements of matrix $A$ are expressed by Equations (6) and (7):

where $S_i$ is the area of tessera $i$ . Equation (6) was defined by Klamt and Schüürmann for a simple sphere partitioned into segments of identical surface area [21]. In C-PCM, the polarization charges are scaled so that the total polarization charge obeys Gauss's law; an effective way to achieve this is to multiply each charge by the factor ${\displaystyle \frac{\epsilon -1}{\epsilon }}$ , yielding new charges given by Equation (8):

Here, $\epsilon$ corresponds to the dielectric constant of the solvent. At the same time, the contribution $\mathsf{\Delta }{G}_{elec}$ is important for polar and charged solutes because of the polarization of the solvent, which is modeled as a uniform medium of a given dielectric constant.

2.2. Free Energy Required for the Construction of a Molecular Cavity in the Solvent

To determine solute–solvent interactions, the molecules are placed in cavities formed by overlapping spheres centered on the solute atoms. Constructing the cavity for the solute in the solvent requires a certain amount of energy, known as the cavitation free energy, a concept introduced by Uhlig [33] and initially applied to determine the solubility of gases in liquids. To determine the cavitation free energy, $\mathsf{\Delta }{G}_{cav}$ , several models have been developed, including the formulation proposed by Reiss and Pierotti [34-36]. This theory describes fluid molecules using hard spheres to evaluate $\mathsf{\Delta }{G}_{cav}$ , which is determined by a polynomial depending on the radius $R_{ms}$ , corresponding to the sum of the radii of the solute and solvent:

The coefficients $K_i$ ( $i=0,1,2,3$ ) are expressed in terms of the properties of the solvent, including its molecular radius and number density, and those of the system (pressure P and temperature T). These can be calculated using Equations (10)–(13):

The term $y={\displaystyle \frac{4\mathsf{\pi }{\mathrm{R}}_s^3{n}_s}{3}}$ is associated with the microscopic terms. The cavitation free energy is a positive term that accounts for the loss of entropy of the solvent due to the reorganization of its molecules in the presence of the solute.

2.3. Free Energies of Dispersion and Repulsion

This section presents the models used to calculate the dispersive and repulsive contributions to the solvation free energy. Using the atom-atom interaction model [37], the dispersion free energy between two molecules A and S is described by Equation (14):

Here, $d_{at}^{\left(n\right)}$ is the dispersion coefficient and $r_{at}$ is the distance between atoms. In real solutions, the description of the relative distances between the atoms of the solute and solvent can be replaced by the distribution function defined in Equation (15):

In Equation (15), $N_t$ is the number of atoms of type $t$ for each solvent molecule, ${\rho }_s$ is the number density of the solvent, and $g_{at}$ is the correlation function, which depends on the interaction distance ( $r_{at}$ ). Using a continuous function allows the summation over S to be replaced by an integral, shown in Equation (16):

Using Equation (16), it is possible to determine the dispersion free energy of a solute in a solvent. In turn, Amovilli [38] proposed a model to calculate the dispersion free energy using a quantum-mechanical treatment:

where $P_A\left(0\mathrm{K}|\overrightarrow{r}\right)$ and ${\omega }_{0\mathrm{K}}^A$ are the densities and energies of transition of a state $K$ of the solute. The term ${\sigma }_B$ is the surface charge density induced by the electric field of the charge distribution $P_A\left(0\mathrm{K}|\overrightarrow{r}\right)$ , which depends on a dielectric constant calculated at imaginary frequencies. Additionally, ${\sigma }_B$ is always smaller than the square of the refractive index ${\eta }_B^2$ , estimated at zero frequency and measured from the electronic transitions in the visible spectrum. However, Mennucci and Amovilli [39] defined these terms using Equation (18):

where ${\Omega }_B={\eta }_B{I}_B$ depends on the ionization potential $I_B$ of the solvent. Here, $E_{0\mathrm{K}}^A\left(\overrightarrow{r}\right)$ is the component of the electric field orthogonal to the surface of the solute cavity at the point $\overrightarrow{r}$ , associated with the transition density $P_A\left(0\mathrm{K}|\overrightarrow{r}\right)$ . Using Equations (18) and (17), the dispersion free energy can be defined as:

The term $V_{0\mathrm{K}}^A$ is the potential associated with the electric field $E_{0\mathrm{K}}^A$ . All terms associated with the dispersion free energy equation are calculated using the restricted Hartree-Fock method.

To calculate the repulsion free energy in the solvation process, Amovilli and Mennucci [39] proposed a model based on the method of Floris [37] and Equation (20):

Here, $A$ and $B$ correspond to the solute and the solvent, respectively. $\overrightarrow{R}$ is an appropriate set of coordinates defining the internal geometry of the complex $AB$ , ${\rho }_B$ is the number density of the solvent, and $g_{AB}$ is a correlation function whose value is zero inside the solute cavity and one outside it. This model is based on the principle that the repulsive force between two interacting molecules follows the Pauli exclusion principle. Based on this principle, the free energy of repulsion increases with the overlap of the two distributions and is related to the density of electrons with the same spin. Thus, the repulsion free energy in the solvation process can be calculated using Equation (21):

Here, $n_{par}^B$ is the number of pairs of valence electrons of the solvent, and $P_{par}\left(\overrightarrow{R}|{\overrightarrow{r}}_2,{\overrightarrow{r}}_1\right)$ and $P_A\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2\right)$ correspond to the adapted pair density matrices. Here, $\overrightarrow{R}$ are the coordinates of the center of the localized orbital containing the pairs of free or bonded electrons of the solute molecule.

2.4. Distribution of a Solute in Two Immiscible Liquids

Consider two immiscible liquids in contact, into which a non-dissociable solute $A$ is incorporated at a concentration within the dilute-solution range [9,40]. The solute $A$ distributes between the two phases to form a solution. If both phases are in equilibrium, the chemical potentials of solute $A$ in both phases are equal and can be written as:

In turn, the chemical potential of solute $A$ in a given phase can be written as:

where $x_{\left(solute,i\right)}$ corresponds to the fraction of solute in equilibrium with one of the phases, and ${\mu }_{\left(solute,i\right)}^0$ is the standard potential of solute $A$ in one of the phases. Using Equations (22) and (23), a relationship can be established to calculate the solute distribution between the two phases:

where the concentration ratio is given by:

Assuming that ${\mu }_{\left(solute,i\right)}^0=\mathsf{\Delta }{G}_{\left(solv,i\right)}^0+\mathsf{\Delta }{G}_{vap\left(A\right)}^0$ , where $\mathsf{\Delta }{G}_{vap\left(A\right)}^0$ is the vaporization free energy of the solute, we can calculate the solvation of $A$ in a given solvent. In this case, the standard chemical potential of solute $A$ in a phase ${\mu }_{\left(solute,i\right)}^0$ is equivalent to the dissolution free energy $\mathsf{\Delta }{G}_{sol}^0$ . Therefore, the distribution constant of the solute in both phases can be written as follows:

Here, $\mathsf{\Delta }{G}_{\left(solv,\beta \right)}^0$ and $\mathsf{\Delta }{G}_{\left(solv,\alpha \right)}^0$ correspond to the solvation free energies of the solute in the two immiscible liquids. Simplifying Equation (26) yields an expression, based on quantum mechanical solvation methods, for calculating the distribution constant $K_{{\displaystyle \frac{o}{w}}}$ of a solute $A$ between two immiscible liquids:

Generally, it is possible to calculate the free energies of phase changes and the solvation free energy. However, determining the free energy associated with the transfer of a solute from one liquid phase to another remains a limitation of quantum mechanical methods used to estimate the chemical potential of a system. Therefore, Equation (26) is particularly suitable for estimating the free energy associated with this phase transfer.

3.1. Quantum Mechanical Calculations

All three-dimensional molecular structures were constructed using the wxMacMolPlt v7.4.5 software [41]. The molecular structures were then optimized in vacuum using density functional theory (DFT), with the PBE-VWN functional and the 6-31G(2p,2d) basis set [42-44].

Using the coordinates obtained from the optimization, solvation calculations in water and 1-octanol were performed using the C-PCM//PBE-VWN/6-31G(2p,2d) method to determine the electrostatic contribution to the solvation free energy. These calculations were based on the theory of Ben-Naim [45], which assumes that a molecule does not change its structural conformation when transferred from vacuum to the solvent medium. Figure 1 shows the solvation process of ethanol in a continuum dielectric medium according to the Ben-Naim principle.

To construct the cavities of the solute molecules, the GEPOL-GB method proposed by Pascual-Ahuir and Tomasi [46,47] was used, together with van der Waals radii for the C, N, H, and O atoms; the radii used are shown in Table 1. For the solvation studies in water and 1-octanol, the dielectric constant, molecular volume, and surface tension were used as input parameters [48,49]. The repulsion and dispersion free energies were calculated using the Amovilli-Mennucci method, which required defining several solvent parameters, including the ionization potential, the square of the refractive index, the density, and the number of valence-shell electrons of the solvent. The ionization potentials of water and 1-octanol were calculated using an SCF//UHF/PBE-VWN optimization procedure with a 6-31G(2p,2d) basis set. All calculations were performed using the GAMESS-US software [50]. The parameters used are shown in Table 2.

3.2. Training Set Selection and Model Validation

The training set of 42 organic molecules used to parameterize the cavitation free energy model (see Table 3) and the set of 50 molecules used for the QSPR models (Models A and B) were selected based on three primary criteria: (i) availability of high-quality experimental solvation free energies and logP values in the literature; (ii) structural and functional diversity, including alcohols, ketones, aldehydes, carboxylic acids, esters, ethers, amines, alkanes (linear and cyclic), and aromatic compounds; and (iii) variation in molecular size and alkyl chain length, to systematically capture hydrophobic contributions. Selecting a diverse training set is essential for developing QSPR models with broad applicability and reliable predictions for unseen compounds.

As an exception, 1,2-ethanediol was included in the training set for the cavitation free energy model. In the logP models, by contrast, alcohol molecules such as 1,2-ethanediol, 1,3-propanediol, and 1,4-butanediol were included in the training set to extend the model's applicability, and glycerol was used to validate this extension.

To prevent overfitting and evaluate the predictive power of the models, two complementary validation strategies were employed. First, leave-one-out cross-validation was performed on the 42 molecules used to parameterize the cavitation free energy model, yielding coefficients of determination of R² = 0.985 for water and R² = 0.982 for 1-octanol, confirming that the models are not overfitted to the training data. Second, external validation was conducted using an independent set of ten molecules (discussed in Section 4.3) not included in the training set, comprising 2-pentanol, n-hexanal, hexylamine, dibutyl ether, propyl acetate, n-butylbenzene, dioxane, anisole, glycerol, and 3-pentanone. The excellent agreement between predicted and experimental logP values for these molecules (discussed in Section 4) provides strong evidence that the models possess genuine predictive capability and are not the result of overfitting. Additionally, all coefficients in the models (see Section 4.1 and Section 4.2) were evaluated for chemical consistency. Positive coefficients were associated with hydrophobic contributions (carbon surfaces and alkyl groups), whereas negative coefficients were associated with polar contributions (oxygen surfaces and hydrogen-bonding groups), confirming that the models capture meaningful structure-property relationships rather than spurious correlations.

4.1. Linear Regression Model to Estimate the Cavitation Free Energy

In general, the solvation free energy is calculated using quantum mechanical methods by means of Equation (1). In this model, the cavitation free energy is not determined accurately by various theoretical models, such as the Claverie-Pierotti method [34-36]. This work therefore develops a model based on a structure–property relationship using the experimental solvation free energies of 42 liquid organic compounds, to allow for more accurate calculation of the cavitation free energy. Table 3 shows the electrostatic, dispersive, and repulsive free energies calculated using the C-PCM and Amovilli-Mennucci methods for this group of 42 organic molecules. Using the Amovilli-Mennucci approximation, the free energy for cavity construction, $\mathsf{\Delta }{G}_{cav}$ , was calculated using Equation (28):

The term $\mathsf{\Delta }{G}_{solv\left(exp\right)}$ corresponds to the experimental solvation free energies reported for these organic molecules in 1-octanol and water. Equation (28) was used by Amovilli and Mennucci to calculate the cavitation free energy: they plotted the free energy associated with cavity construction as a function of the cavity areas ( $A_{cavity}$ ) of the molecules, obtaining a linear relationship of the form:

In Equation (29), the terms $\alpha$ and $\beta$ are parameters obtained by linear regression. This model has been proposed in other solvation studies [51,52]. Here, we propose an equation to calculate the free energy for cavity construction based on molecular descriptors calculated with the C-PCM and GEPOL-GB methods included in the GAMESS-US package. The linear regression model is shown in Equation (30):

where $n=4$ and $m=9$ . The term $i$ corresponds to the different atoms in the solute molecule ( $i=$ C, H, N, and O). The factor ${\sigma }_i$ is associated with each atom, and $k_j$ is a factor related to each functional group present in the molecule. This equation, obtained by linear regression, includes the surface areas of the atoms in the solute molecule and a factor associated with the functional groups present in the molecule. The term $X$ is equal to one if the functional group is present and zero otherwise.

The parameter $b$ is a constant term obtained by linear regression, and $A_i$ is the surface area of atom i obtained using the GEPOL-GB method [46]. Models similar to Equations (29) and (30), based on statistical mechanics, have been reported by Marenich, Cramer, Truhlar, and Kollman [53-55].

The parameter values obtained by linear regression to estimate the free energy required for cavity construction are given in Table 4. This model was parameterized using the root mean square error $\left(RMSE\right)$ , defined by:

The parameterization was performed for water and 1-octanol, with coefficients of determination of $R^2=0.990$ and $R^2=0.986$ , respectively. In the case of 1-octanol, the coefficient $k_j$ was not determined for the acid functional group.

As presented in Table 5, the solvation free energies in water and 1-octanol were obtained using the free energy of cavity formation, as estimated by Equation (30). The results are in good agreement for the 42 molecules included in the parameterization process.

As illustrated in Figure 2 and Figure 3, the estimated solvation free energy values are plotted against the experimental solvation free energy in water and 1-octanol. For the multivariate model, the linear regression curve demonstrates a strong correlation for the solvation free energies in both water and 1-octanol. The errors are generally larger for 1-octanol solvation energies than for water solvation energies, as evident from the RMSE values in Table 4: RMSE = 1.25 kJ mol⁻¹ for water and RMSE = 0.67 kJ mol⁻¹ for 1-octanol. Although the RMSE is smaller for 1-octanol, the relative errors (see values in Table 5) show that the deviations for 1-octanol are more systematic and depend on the functional group. This indicates that the use of macroscopic solvent parameters for 1-octanol introduces systematic errors, because this solvent exhibits a complex molecular structure that is not adequately represented by a uniform dielectric continuum. 1-Octanol forms dynamic aggregates through hydrogen bonding and hydrophobic interactions, creating a heterogeneous environment that affects solute partitioning in ways that simple continuum models cannot capture. This is particularly evident for molecules with polar functional groups (e.g., alcohols, ketones), where specific interactions with the 1-octanol hydroxyl group are important.

The lowest correlation was obtained for the solvation free energy in water. This is due to the 6-31G(2p,2d) basis set used to estimate the electrostatic contribution and the minimal basis set used in the Amovilli-Mennucci method to determine the dispersion and repulsion free energies. In other studies, these solvation free energies have been estimated using additional molecular descriptors, such as hydrogen-bond acceptor and donor indices, which has improved the estimation of these properties [56,57].

4.2. Prediction of the Partition Coefficient K o / w

Using the estimated solvation free energies and Equation (27), the partition coefficient was determined. Table 6 shows the estimated partition coefficient values, which vary considerably from the reported experimental values.

The coefficient of determination obtained was $R^2=0.974$ , using the solvation free energy values calculated in this work (see correlation in Figure 4). These results show a better correlation than those reported by Amovilli-Mennucci using their model to estimate the cavitation free energy.

At the same time, Table 6 presents the partition coefficients with their respective relative errors, showing that the greatest estimation error occurred for 2-propanol and molecules containing the carbonyl functional group. Additionally, the amines exhibit high relative errors in the estimated partition coefficients. Nevertheless, the results obtained for the partition coefficient using solvation free energies estimated with the cavitation free energy model are better than the values obtained by Amovilli-Mennucci. It is also important to note that the distribution of a solute between two immiscible solvents is fundamentally a transfer process; therefore, using solvation free energies is not entirely suitable for estimating partition coefficients, a limitation that should be addressed using a model that determines the transfer free energy directly. The 6-31G(2p,2d) basis set used to calculate the electrostatic contribution to the solvation free energy remains appropriate; however, a systematic source of error is the minimal basis set used by the Amovilli-Mennucci method to calculate the dispersion and repulsion free energies.

As shown in Table 6, molecules with significant van der Waals contributions to their solvation in 1-octanol, such as 2-propanol and dimethyl ether, exhibit the largest deviations. In these cases, the minimal STO-3G basis set is too limited to describe the anisotropic electron density distributions necessary for accurate dispersion interactions. This limitation is particularly pronounced for molecules containing secondary alcohol groups, in which the steric environment surrounding the hydroxyl group modulates dispersion interactions with the solvent in a manner that requires a more flexible basis set. In future work, basis sets of at least triple-zeta quality with polarization functions, such as 6-311+G(d,p), will be employed to improve the accuracy of these non-electrostatic contributions. It is important to note that the deviations observed in Table 6 are not random but reflect a fundamental limitation of the approach based on the difference in solvation energies (see Equation (27)). This method assumes that the transfer free energy between phases is exactly equal to the difference in solvation energies, a simplification that ignores activity effects, conformational changes, and molecular association in the organic phase. Additionally, although the GEPOL-GB method for cavity construction is computationally efficient, it relies on fixed van der Waals radii and a single molecular conformation optimized in vacuum. Consequently, this approach does not account for conformational flexibility, which is particularly problematic for molecules capable of adopting multiple low-energy conformations with different solvent-accessible surface areas. For example, 2-propanol can rotate around its C–O bond, altering the exposure of its hydroxyl group to the solvent. 1,2-Ethanediol, by contrast, shows an excellent prediction because its conformation is relatively constrained by intramolecular hydrogen bonding, whereas flexible molecules like 2-propanol and methyl butanoate can adopt multiple conformations, and using a single optimized structure fails to capture the conformational ensemble that contributes to solvation.

In this paper, we also propose two quantitative structure–property relationship (QSPR) models for estimating the partition coefficient using properties obtained with the C-PCM method. These QSPR models incorporate the electrostatic free energy calculated in 1-octanol and water, the atomic surfaces of the atoms composing the molecules (determined using the GEPOL-GB method), and factors associated with the functional groups present in the molecule. In this case, 50 molecules were used for parameterization. The first QSPR model, defined as Model A, is shown in Equation (32):

where $n=4$ and $m=9$ . The term $i$ corresponds to the different atom types in the solute molecule ( $i=$ C, H, N and O). The factor ${\sigma }_i$ is the atomic tensor associated with each atom and $k_j$ is a factor related to each functional group present in the molecule. The terms $\mathsf{\Delta }{G}_{elec\left(water\right)}$ and $\mathsf{\Delta }{G}_{elec\left(1-octanol\right)}$ correspond to the electrostatic free energies of solvation and the term $X$ is equivalent to one if the functional group is present and zero otherwise.

Model A achieves a correlation coefficient of $R^2$ = 0.991, making it simpler than the model of Equation (30), which requires greater computational cost for estimating the free energy of dispersion and repulsion. This model is also more accessible and better suited to large molecules.

All parameters of Equation (32) were obtained using the root mean square error $\left(RMSE\right)$ , defined by Equation (33):

Table 7 contains the values obtained by multivariate regression of the coefficients associated with Model A. Unlike the model that uses the cavitation free energy, Model A directly yields the distribution coefficient. The RMSE was 0.16, and the statistical analysis shows that the most influential factors are the presence of aromatic and amine functional groups. Model A is parameterized directly using experimental logP values, which incorporates the transfer free energy between phases and avoids the systematic errors inherent in the solvation energy difference approach. The chemical significance of the functional group coefficients in Model A reflects the distinct solvation environments in water versus 1-octanol: negative coefficients for aromatic groups (−3.168) and amines (−3.789) indicate their strong preference for the organic phase due to π-π stacking and hydrophobic effects, while the positive coefficient for acids (0.189) reflects their greater affinity for water through hydrogen bonding. The ketone coefficient (−1.457) reflects the balance between carbonyl hydrogen bonding with water and dispersion interactions with the alkyl chain of 1-octanol, whereas the ether (−1.208) and ester (−0.783) coefficients reflect their intermediate polarity and their ability to act as hydrogen bond acceptors in both phases. The atomic surface coefficients enhance the performance of Model A by accounting for molecular size and shape effects, with the positive carbon surface coefficient (0.086) capturing hydrophobic contributions and the negative oxygen surface coefficient (−0.063) reflecting polar solvation effects. This combination of direct transfer free energy modeling, empirically calibrated functional group contributions, and atomic surface descriptors allows Model A to achieve superior correlation (R² = 0.991, RMSE = 0.16) while maintaining computational simplicity and chemical interpretability. Table 8 shows the estimated partition coefficient values and their absolute relative errors using Model A, where dimethyl ether and 2-propanol show the greatest relative error. Figure 5 shows the estimated partition coefficient values for Model A relative to the experimental partition coefficient.

Finally, we present a second model to determine the partition coefficient of the 50 organic molecules, defined as Model B. Model B uses properties such as the non-bonding free energy, given by Equation (34):

This multivariate regression model includes the atomic surfaces $A_i$ of the atoms composing the molecules and the coefficients associated with their functional groups. Model B is shown in Equation (35):

Similarly, the coefficients in Equation (35) were obtained using the experimental partition coefficient values in the multivariate regression. Figure 6 shows the estimated values for Model B versus the experimental partition coefficient values. Model B achieves a correlation coefficient of 0.986 and a Root Mean Square Error (RMSE) of 0.17 relative to the experimental data. Table 9 shows the parameter values obtained for Model B. In this model, no coefficients are associated with the aldehyde, alcohol, ester, and cycloalkane functional groups, which makes it a simpler model, although determining coefficients for these functional groups would have improved the multivariate correlation. However, Model B has a higher computational cost because it involves the calculation of dispersive and repulsive free energies. The systematic deviations observed for alcohols (2-propanol, relative error 6.14), amines (diethylamine, relative error 0.55), and carbonyl-containing compounds (acetone, relative error 0.84; methyl acetate, relative error 1.64) indicate that hydrogen bonding and specific intermolecular interactions are not fully captured by the current models, reflecting the inherent limitations of continuum solvation methods that represent solvents as uniform dielectric media. Including additional descriptors, such as hydrogen-bond donor/acceptor ability, molecular dipole moment, or polar surface area, would likely improve predictive capability by explicitly encoding information about the solute's capacity for specific interactions. While the atomic surface coefficients in these models partially capture polar solvation effects, they do not distinguish between hydrogen bond donors and acceptors, nor do they account for interaction strength or directionality.

Additionally, Table 10 presents the logP values predicted by Model B for 50 organic molecules; the largest relative errors were observed for molecules containing ester, hydroxyl, and aldehyde functional groups, as the corresponding parameters were not included in the analysis. Furthermore, both QSPR models are limited to molecules containing a single functional group within their structures.

The linear regression models assume additive contributions from atomic surfaces and functional groups. However, the results indicate that this additivity assumption breaks down for molecules with intermediate alkyl chain lengths (e.g., 2-butanone), where the balance between hydrophilic and hydrophobic contributions is most sensitive to the local molecular environment. This suggests that non-linear effects, such as synergistic interactions between the carbonyl group and the alkyl chain, are not fully captured by the current model. Future studies should explore non-linear regression techniques, such as random forests or neural networks, using a larger and more diverse training set to better capture these complex structure-property relationships. Furthermore, $K_{o/w}$ describes the equilibrium distribution of a solute between two immiscible liquid phases and is fundamentally a property of the transfer process between solvents, not a property of solvation from vacuum. While solvation free energies (ΔGsolv) describe the transfer from vacuum to a single solvent, the transfer free energy (ΔGtransfer) describes the transfer between solvents and is directly related to logP. Using solvation free energies to estimate logP (Equation (27)) introduces systematic errors because it assumes that the difference in solvation free energies (both referenced to vacuum) equals the transfer free energy between phases. This assumption neglects the differences in activity coefficients between solvents, conformational changes induced by different solvent environments, and specific solute-solvent and solvent-solvent interactions that are not captured by continuum models. These limitations explain why the solvation energy difference approach in Table 6 shows larger deviations (R² = 0.974) compared with Model A (R² = 0.991), which is parameterized directly against the experimental logP values and effectively models the transfer free energy.

Additionally, the dispersive and repulsive free energies in Model B account for van der Waals interactions and Pauli exclusion effects that are particularly significant in 1-octanol, where the long alkyl chain provides extensive dispersion interactions with hydrophobic regions of the solute, while repulsive terms capture excluded-volume effects governing molecular packing and cavity formation in both phases. However, these contributions are calculated using a minimal basis set (STO-3G) in the Amovilli-Mennucci method, which systematically underestimates dispersion interactions for molecules with diffuse electron densities and fails to capture effects that are crucial for accurate logP prediction, limiting their overall contribution to predictive accuracy. Including these terms in Model B yields only a marginal improvement in correlation (R² = 0.986) compared with Model A (R² = 0.991), while introducing additional computational cost and potential error propagation—demonstrating that the dominant factors for logP prediction are electrostatic interactions and molecular surface descriptors, rather than explicit dispersion-repulsion calculations. The systematic errors in the dispersion-repulsion terms also explain why the solvation energy difference approach (see values in Table 6) exhibits larger deviations for molecules in which van der Waals interactions dominate solvation in 1-octanol, such as alkanes and ethers, whereas Model A effectively circumvents this limitation through direct parameterization against experimental data. Consequently, despite their theoretical importance, the practical contribution of dispersive and repulsive free energies to predictive accuracy is limited by current methodological constraints, making Model A the most computationally efficient and accurate choice for routine logP prediction.

4.3. Validation of QSPR Models with Additional Molecules

The results presented in Table 11 demonstrate the predictive performance of Models A and B for ten additional organic molecules that were not included in the training set. These molecules encompass a chemically diverse range of functional groups, including secondary alcohols (2-pentanol), aldehydes (n-hexanal), amines (hexylamine), ethers (dibutyl ether, dioxane, and anisole), esters (propyl acetate), aromatics (n-butylbenzene), and polyols (glycerol). Both models show generally good agreement with the experimental logP values, with Model A achieving a correlation coefficient of R² = 0.991 and Model B achieving R² = 0.986, confirming the robustness and general applicability of the QSPR approach across different chemical families. For most compounds, the predicted logP values are within 0.2–0.3 log units of the experimental values, with n-butylbenzene showing nearly perfect prediction by both models (Model A: 4.26, Model B: 4.38, Exp: 4.38) and glycerol showing the largest deviations (Model A: −2.02, Model B: −2.59, Exp: −1.76), attributable to the strong hydrogen bonding capabilities of the triol, which are not fully captured by continuum solvation models.

The prediction for hexylamine by Model A (4.98 versus 2.06 experimental) reveals a significant overestimation, suggesting that the amine functional group coefficient in Model A may require further refinement for primary amines with longer alkyl chains, whereas Model B provides a much more accurate prediction (2.08 versus 2.06), demonstrating the value of including dispersive and repulsive contributions for molecules where van der Waals interactions dominate solvation in the organic phase. For ethers and esters, both models perform well. Dibutyl ether was predicted to have logP values of 3.27 (Model A) and 4.01 (Model B), compared with the experimental value of 3.21, whereas propyl acetate was predicted to have logP values of 1.12 (Model A) and 1.14 (Model B), compared with the experimental value of 1.24. These results indicate that the atomic surface descriptors effectively capture the hydrophobic contributions of alkyl chains. The prediction for anisole (Model A: 2.03, Model B: 1.88, Exp: 2.11) confirms that the aromatic coefficient (−3.168 in Model A) appropriately accounts for the π-π interactions and hydrophobic character of the methoxybenzene structure. For dioxane, both models correctly predict a negative logP value (Model A: −0.14, Model B: −0.67, Exp: −0.27), reflecting the hydrophilic nature of the cyclic ether due to its two oxygen atoms, although Model A provides better accuracy owing to its direct parameterization against experimental data. The 3-pentanone prediction (Model A: 1.96, Model B: 2.06, Exp: 2.03) demonstrates that the ketone functional group coefficient successfully captures the balance between carbonyl hydrogen bonding and hydrophobic alkyl chain effects for symmetrical ketones. These results collectively confirm that Model A, despite its computational simplicity, provides predictions comparable to or better than Model B for most compound classes, with the exception of long-chain amines, where dispersion contributions become important—highlighting the complementary nature of both approaches across different chemical spaces.

4.4. Deviations and Limitations

The deviations observed in Table 6 between the estimated and experimental logP values are systematic rather than random and can be attributed to several factors. First, it is important to note a fundamental limitation of the solvation energy difference approach illustrated in Equation (27): this approach assumes that the free energy of transfer is exactly equal to the difference in solvation energies, an assumption that fails to account for activity effects, conformational changes, and molecular association in the organic phase. This limitation is most severe for molecules such as 2-propanol and dimethyl ether, which can adopt different conformations or interact specifically with the solvent. Second, the use of a minimal basis set for the dispersion and repulsion terms is a significant concern: the Amovilli-Mennucci method uses a minimal basis set (STO-3G) by default, which inadequately describes the anisotropic electron density distributions necessary for accurate dispersion interactions. This is particularly relevant for molecules in which van der Waals interactions predominate over other forces during solvation in 1-octanol. Third, the GEPOL-GB cavity model, which uses fixed van der Waals radii and a single vacuum-optimized conformation, does not account for conformational flexibility or averaging, presenting a significant challenge for molecules that exhibit multiple low-energy conformations. Fourth, the assumptions underlying linear regression models are relevant: the additivity assumption becomes invalid for molecules with intermediate alkyl chain lengths, owing to the sensitivity of their hydrophilic and hydrophobic contributions to the local molecular environment. Furthermore, macroscopic solvent parameters for 1-octanol fail to capture its microheterogeneous structure (monomers, dimers, aggregates), leading to systematic errors for molecules with polar functional groups.

Finally, it is important to emphasize that Model A (Equation (32)) substantially mitigates these limitations by being parameterized directly against experimental logP values, achieving a considerably better correlation (R² = 0.991, RMSE = 0.16) than the solvation energy difference approach (R² = 0.974). By contrast, the Amovilli-Mennucci method systematically underestimates dispersion interactions for molecules with diffuse electron densities and fails to capture effects that are crucial for accurate logP prediction, limiting its overall contribution to predictive accuracy; including these terms in Model B yields only a marginal improvement, with a correlation coefficient of R² = 0.986. Additionally, a limitation of the current QSPR models is their use of binary indicators for functional groups, whereby a value of 1 is assigned if the functional group is present and 0 otherwise, regardless of the number of identical groups in the molecule. For example, 1,2-ethanediol receives the same alcohol indicator as methanol or ethanol, despite having two hydroxyl groups that substantially increase its hydrogen bonding capacity and hydrophilicity. This simplification is particularly problematic for polyols such as glycerol, which contains three hydroxyl groups and shows larger deviations (Model A: −2.02, Model B: −2.59, Exp: −1.76) compared with monohydric alcohols, suggesting that the models underestimate the hydrophilicity of polyfunctional compounds. Using functional group counts (e.g., number of OH groups, NH groups, or carbonyl groups) rather than binary presence would allow the models to capture the additive contributions of multiple identical groups to overall partitioning behavior, improving predictions for polyfunctional molecules. Furthermore, for compounds containing multiple different functional groups simultaneously (e.g., amino alcohols, hydroxy acids, keto-esters), interaction terms between functional groups could be incorporated to account for non-additive effects, such as intramolecular hydrogen bonding, which can significantly modify solvation properties compared with the sum of individual group contributions. The binary indicator approach was chosen in this work for its simplicity and interpretability, as well as the limited number of polyfunctional compounds in the training set.

In this work, the partition coefficients of several organic molecules were estimated using quantitative structure–property relationship (QSPR) models. Constructing these models required only the electrostatic free energy obtained by the C-PCM//PBE-VWN/6-31G(2p,2d) method and the dispersive and repulsive free energies calculated by the Amovilli-Mennucci method. Of the three proposed models containing parameters calculated with the C-PCM method, Model A achieved the highest coefficient of determination (R² = 0.991), with an RMSE of 0.16. This model relates the electrostatic free energy and the terms associated with the functional groups present in the molecule to the partition coefficient. Additionally, the QSPR models obtained in this work to estimate the free energy for cavity construction in water and 1-octanol were acceptable, with correlation coefficients of R² = 0.990 and R² = 0.986 for 1-octanol and water, respectively.

While a direct comparison with alternative solvation models such as SMD, COSMO-RS, or IEFPCM would further validate the robustness of these conclusions, the primary objective of this work is to demonstrate that C-PCM-derived descriptors, when combined with QSPR modeling, can achieve high predictive accuracy for logP (R² = 0.991) as a revision and extension of the original Amovilli-Mennucci methodology. The C-PCM//PBE-VWN/6-31G(2p,2d) approach was selected because it provides a favorable balance between accuracy and computational cost. Moreover, its wide availability in the GAMESS-US package ensures accessibility for the broader research community. Although SMD and COSMO-RS incorporate more refined parameterization of non-electrostatic contributions, the excellent correlation achieved by our model suggests that the electrostatic and surface descriptors derived from C-PCM capture the essential physics governing the partitioning process. Future work will explore the application of SMD, COSMO-RS, and IEFPCM to further improve the predictive capability of these QSPR models and evaluate their general applicability across different solvation methodologies. Finally, the predictive accuracy of the proposed models would be substantially enhanced by expanding the training dataset to include a wider range of chemical functionalities, particularly molecules with multiple functional groups, heterocycles, and halogenated compounds. Further external validation using larger independent datasets not included in the parameterization process would strengthen the reliability and generalizability of the models.

The author is solely responsible for conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing—original draft preparation, writing—review and editing, and visualization of the article. The author has read and agreed to the final version of the manuscript.

The data supporting the findings of this study are available within the article and from the author upon reasonable request.

The author declares no conflicts of interest.

The study did not receive any external funding and was conducted using only institutional resources.

The author gratefully acknowledges the Universidad de Carabobo, Facultad Experimental de Ciencias y Tecnología, and the Laboratorio de Simulación Molecular para la Química Experimental (SIMOLQUIMEX) for providing the computational resources and institutional support necessary to carry out this work.

In accordance with COPE guidance on "Authorship and AI Tools" (2023), DeepSeek (DeepSeek-V3/DeepSeek-R1) was used solely to edit existing text for language polishing, grammatical correction, stylistic refinement, and improved clarity, readability, and coherence. It was not used as an author, nor was it used to generate scientific content, data, methodology, interpretations, or conclusions. All AI-assisted edits were performed under my supervision and were critically reviewed, verified, and approved before incorporation. As the sole author, I take full responsibility for the content, accuracy, integrity, and originality of the work, including all data, calculations, interpretations, and conclusions.