Free energy calculations and solubility in water of organic molecules: a numerical relation through molecular dynamics

ABSTRACT Solubility in water, free energy of solvation (ΔG ) and free energy of hydration (ΔG ) were determined for 86 organic liquids using Molecular Dynamics simulations in order to find out the relation among these properties. Both, correlation between ΔG and solubility and correlation between ΔG and solubility are not clearly defined. However, by defining =ΔG G , it is found a numerical relation between and solubility. This result was verified for pentane, ethanol and, hexanol of three different force fields (TraPPE-UA, CHARMM and, GROMOS). Also, it was studied the relation of these properties by increasing hydrophobicity in a functional group through the number of carbon atoms of the hydrocarbon chain. In this last case, the mentioned relation was also found.


Introduction
Solubility is a property which determines how much a substance can stay in a solution without precipitation and, therefore, is defined as the maximum amount of a solute that can be dissolved in a solvent under given physical conditions (pressure, temperature, pH, etc.). The solubility of a compound can change when pressure, temperature, and/or compositions change [1]. In the case of water, the more polar the compound is the more soluble it is. The solubility of a constituent within a multi-component mixture may be orders of magnitude lower than the aqueous solubility of the pure constituent in water. For many compounds, the solubility increases with temperature, a compound may have considerably different solubility depending on the investigated solvent system and additionally, the solubility of different compounds may be totally different in a specific solvent. In the case of oils, their solubility in water varies between them because they are a complex mixture of compounds [2].
The solvent is, in general, a liquid which can be a pure substance or a mixture of two liquids. It is possible also speak of solid solution, but rarely of solution in a gas. The extent of solubility goes from infinitely soluble (fully miscible) such as ethanol in water to poorly soluble, such as silver chloride in water. The term insoluble is applied to poorly or very poorly soluble compounds. Solubility occurs under dynamic equilibrium, which means that solubility results from the simultaneous and opposing processes of dissolution and phase joining. Solubility equilibrium occurs when the two processes proceed at a constant rate. Under certain conditions, equilibrium solubility may be exceeded to give a so-called supersaturated solution, which is metastable [3].
Solubility is used to take very important decisions from the earliest stages of drug discovery, throughout the entire process of development and up to formulation. For many products, crystallization is used for purification as well as particle formation. In the crystallization of active ingredients, the solubility curve helps to choose a suitable crystallization process and determines the yield. Thus, knowledge of the solubility is essential for the design of the crystallization process. A widely accepted and accurate method for measuring the solubility is through equilibration of a suspension, followed by an assessment of the solution composition, from which the solution concentration can be determined [4][5][6]. Experimentally the preferred method for determining the solubility of a solute at various temperatures has been to find the temperature at which a solution of known composition is saturated, that is, at which the solute begins to crystallize out of solution [7].
The solubility theories are based on the condition that, for the occurrence of solubility of two liquids, it is necessary that the intermolecular interactions between the molecules of one component (A-A) and the molecules of a second component (B-B) are of the same order of magnitude and, therefore, can be broken to form A-B interactions [8]. Hildebrand established the first theory of solubility parameters and described it as being the square root of the cohesive energy density [9]. In the same way, for polymer dissolution, it is necessary that the solubility parameters of the solvent and solute are close to each other. However, this solubility parameter does not discriminate the kind of intermolecular interactions and in some cases predicted incorrect results [10].
To improve the Hildebrand's methodology, Hansen proposed a new treatment based on the three different types of intermolecular interactions: dispersive interactions (nonpolar interactions, δD), polar interactions (between permanent dipoles, δP), and hydrogen bonds (δH). Therefore cohesive energy density can be described as the square root of the sum of the quadratic of the three different kinds of interactions. Thus, for Hansen, a solvent could be represented as a point in a three-dimensional space, with coordinates δD, δP, and δH. Similarly, a solvent mixture can also be represented as a point in three-dimensional space, but in this case with coordinates d m D , d m P , d m H [10]. The theoretical bases for the prediction of solubilities are provided by thermodynamics and statistical mechanics. In many cases, statistical mechanics yields the chemical potential as an equation including several terms. While one or more terms may be fully determined theoretically, the parameters in other terms must be evaluated empirically at present. Then the solubility is obtained by equating the chemical potential for each species in different phases or by minimising the total Gibbs energy at constant pressure and temperature [11].
For the prediction of gas solubility, molecular simulation has proven successful in both organic solvents and water. The quantitative accuracy of these predictions crucially depends on the quality of molecular models for solute and solvent species. Through calibration against experimental data, solute models can be optimised to reproduce gas solubility in a particular solvent. Also, ensuring a correct representation of solubility for each pure solvent amounts to a considerable effort in model development and validation, particularly for the simulation of mixed solvents with more than two components [12].
By molecular simulations, recently, it has been presented a methodology to determine the solubility of crystal organic molecules in solution based on the equilibrium free energy difference between the solvated solute and its crystallized state at the crystal surface kink site [13].
Recently, a review of computational methods to determine the solubility of pharmaceutical molecules was published. Though the main focus is on classical molecular dynamics simulations, the authors also present other computational techniques, such as the COSMO method, and also discuss Flory-Huggins theory and solubility parameters This review gives insights into the methods for calculating solubility instead experiments [14].
Also, in recent years solubility has been determined through density profiles of equilibrated molecular mixtures in order to reparameterise and improve force fields by making them capable to reproduce not only dielectric constant, surface tension, density and, solubility, but also, the transition between molecules with the same functional group soluble and insoluble in water (e.g. propanol and butanol) [15,16].
A hybrid approach for solubility prediction of gases in mixed solvents which combines experimental data for Henry's law constants in pure solvents with molecular simulation in order to estimate non-ideal contributions to solubility in mixtures has been also proposed, however, in this method, the solute model is not optimised so that methodology only describes in good agreement with experimental data the solubility [12].
An attempt to predict the mutual solubility of water and dipropylene glycol dimethyl ether by directly calculating the Gibbs free energy change of mixing was reported. In that work authors mentioned that the challenges are mainly in force field quality and simulation efficiency, they modified the combining rules to effectively include the polarisation effect and extended the simulations to the limit that their computational resources could provide. The modification of force field was validated by calculating hydration free energies of dimethyl ether and diethyl ether. The average uncertainties in predicted free energy changes are less than 0.8 kJ/mol. Unfortunately, this unprecedented low uncertainty was still too large in comparison with the absolute values of ΔG mix [17].
The solvation free energy gives the free energy change associated with the transfer of a molecule between ideal gas and solvent at a certain temperature and pressure. While solvation free energies (ΔG solv ) in general and hydration free energies (ΔG hyd , solvation in water) in particular might not seem to have far-reaching implications, in fact, research can benefit from their prediction because such solvation free energies are related to a broad range of physical properties such as infinite-dilution activity coefficients, Henry's law constants, solubilities, and distributions of chemical species between immiscible solvents or different phases [18].
Solvation free energies not only tell us how much a molecule prefers one phase over another but also can provide insight into how a solvent behaves in different environments. For example, water solvates molecules of opposite polarity differently because of its inherent asymmetry [19], surfaces also have asymmetric effects on ion pairing that depend on the curvature of the surface [20], and the molecular geometry and chemical environment affect hydrophobic solvation [21]. Although ΔG solv and ΔG hyd can be difficult to measure experimentally, they can be calculated to a precision of better than 0.4 kJ/mol, even with a relatively modest investment of simulation time, for relatively diverse small neutral molecules such as those seen in the FreeSolv database [22] of hydration free energies and in recent blind challenges such as the Statistical Assessment of the Modeling of Protein and Ligands challenges. These challenges aim to improve the quality of predictive computational tools in drug design and have leveraged solvation free energies to help drive improvements in modelling.
In the present work, it is shown a clear evidence of the relationship among solubility in water, ΔG solv and ΔG hyd , where ΔG solv refers to free energy calculated for a molecule A solvated by a system of A molecules and ΔG hyd refers to free energy calculated for a molecule A solvated by water. Solubility was determined through density profiles of equilibrated water mixtures and free energy calculations were performed by thermodynamic integration.

Free energy calculation
The Gibbs free energy is defined as: where A is the Helmholtz free energy, p is the pressure, and V is the volume. Obtaining Gibbs free energy difference of a state X to a state Y is performed by considering both G and A dependent on λ, this parameter indicates how the system smoothly evolves from an X state to a Y state, namely, the X state is represented by a λ= 0 and the Y state by λ = 1, the intermediate states are meaningless but turn on smoothly the Y state. The explicit dependence of this parameter is through the force field; Lennard-Jones, Coulomb, bonds, and angles potentials, constraints, dihedrals, kinetic energy, etc. So, to obtain the values from the simulation of G and A is by integrating, over the state X to Y or explicitly from 0 to 1 in the λ dependence, the derivatives respect to λ of both energies [23].
In practice, to evaluate the difference in the free energy between subsequent states it was calculated the derivative of the total energy (the so-called Hamiltonian, including both potential and kinetic) in the system. This methodology is known as Thermodynamic integration (TI) and is one of the most widely used free energy calculation approaches [24]. This method covers a transition between an initial and a final state, and similarly to Free Energy Perturbation calculations uses a coupling parameter λ to describe the transition between these two points [25,26]. By introducing several intermediate λ values between 0 and 1, the energy function as a function of λ would be defined as in The potential energy U(λ) in each λ-state can be calculated as an ensemble average over configurations sampled from an MD simulation. The change in free energy between the initial and final states 0 and 1 can then be computed from the integral of the ensembleaveraged derivatives of potential energy over the coupling parameter λ, where ΔG is the free energy and U(λ) is the potential energy expressed as a function of the coupling parameter λ [27].

Computational details
Gibbs free energy of solvation DG solv and hydration DG hyd were calculated in an NPT ensemble, where N is the molecule number, P is the pressure, and T the temperature. The systems to obtain both energies were built with 1 molecule of the compound to be studied solvated with 499 molecules of the same compound or hydrated with 499 water molecules to obtain DG solv or DG hyd , respectively. One molecule solvated by this number of molecules is enough to assume infinite dilution and perform a free energy calculation properly [28].
To keep constant temperature and pressure, Nose-Hoover and Parrinello-Rahman were applied with relaxation times of t t =0.2 ps and t p =5.0 ps, respectively. A cut-off radius of 1 nm and 10 ns of simulation for equilibrating system and the last 4 ns were used for calculations.

Solubility
The solubility of a solvent i in a solvent j was obtained as, where r i and r j are the density profile for the solute and solvent, respectively. M i is the solute molar mass (kg/mol). The calculation of solubility was done by considering the density rate of the solute and solvent plateau multiplied by 1000 g/L [29]. Through equilibrated density profiles it is possible to determine the solubility by analysing the mass amounts mixed in equilibrium. In Figure 1(A) it can be observed a homogeneous mixture where profiles show that both, water (red line) and the set of molecules represented by the black line are present along the Z direction with the same density; this result denotes a totally miscible mixture. In Figure 1(B) it can be seen that molecules represented by the black line show certain solubility in water, namely, in the region where water is predominant (≈6 nm to ≈10 nm) molecules are present with a density of 50 ± 5 kg/m 3 , since the water density showed in the profile is 950 ± 5 kg/m 3 the solubility is 52.6 ± 5.05 g/L and it will be referred as partial. Finally, Figure 1(C) shows a totally immiscible mixture without any mixture in both regions, in this case, solubility is zero and it will be referred to as insoluble. It is important to mention that density profiles are an average of many configurations and in some cases, for systems reported as insoluble there are a few molecules solvated in the solvent, however, this number of molecules is not significant to consider the system as partially soluble. Solubilities for the 86 organic liquids were determined by density profiles in such a way.

Computational details
To obtain solubility an NP N AT semi-isotropic ensemble was established [30], where N is the particle number, P N is the normal pressure (z-direction), A is the interfacial area, and T is the temperature. The temperature and normal pressure were set at 300 K and 1 bar, respectively. The Nose-Hoover and Parrinelo-Rahman with relaxation times of t t =0.5 ps and t p =1.0 ps kept constant temperature and normal pressure, respectively. An interfacial area of ≈ 15 nm 2 and a truncation distance of 1.5 nm were set for all simulations. For this calculation, it is established a liquid/liquid system, with the interface along the Z direction. As initial configuration, the molecular systems are not mixed. For these calculations it is necessary to have two liquid systems which allow us to obtain density profiles well defined in the plateau in order to determine density; in the case of water 2500 molecules is enough, for the organic molecules it depends on the molecule size (i.e. N=369 for 2- hexanone and N=851 for nitromethane), since X = Y ≈ 3.9 nm, Z dimension for each system has to be at least equal to X. Simulation time was 100 ns and the last 10 ns were taken for results production [15,16,31].
All simulations were performed using version 5.0.7 of GROMACS [32] where the equations of motion are solved using the leap-frog algorithm with a time step of 2 fs. Periodic boundary conditions are used in all directions. The Particle-Mesh Ewald (PME) and the linear constraint solver algorithm (LINCS) were established [33].

Force fields
This section describes the interaction potentials and the parameters used to perform the molecular simulations, as well as the simulation details.
2.3.1. Optimised potentials for liquid simulations force field (OPLS) [34][35][36] The functional form is based on the OPLS/AA force field, where the intermolecular potential contains Lennard-Jones and Coulomb contributions.
where r ij is the separation distance between two particles, s ij is the distance at which U LJ (r)=0, e ij is the well depth. θ is the angle distance and ϕ is the dihedral angle; k r and k u are the constants for the harmonic oscillators and V 1 , V 2 , V 3 and V 4 are the constants for the torsional function.

Chemistry at HARvard macromolecular mechanics (CHARMM) [37,38]
The functional form CHARMM where r ij is the separation distance between two particles, σ is the distance at which U LJ (r)=0, e is the well depth. The energy function accounts for the bond stretches where k b is the bond force constant and b − b 0 is the distance from equilibrium that the atom has moved, the bond angles where k u is the angle force constant and u − u 0 is the angle from the equilibrium between 3 bonded atoms, the dihedrals (torsion angles) where k f is the dihedral force constant, n is the multiplicity of the function, ϕ is the dihedral angle and δ is the phase shift. the impropers, that is out of plane bending, where k v is the force constant and v − v 0 is the out of plane angle. The Urey-Bradley component (cross-term accounting for angle bending using 1,3 nonbonded interactions) comprises the term, where k f is the respective force constant and ϕ is the distance between the 1,3 atoms in the harmonic potential. Nonbonded interactions between pairs of atoms (i, j) are represented by the last two terms. By definition, the nonbonded forces are only applied to atom pairs separated by at least three bonds.

Transferable potentials for phase equilibria (TraPPE) [39,40]
The Transferable Potentials for Phase Equilibria force fields is a collection of functional forms and interaction parameters useful for modelling complex chemical systems, which functional form is: where r ij , e ij , s ij , q i , and q j are the separation, LJ potential well depth, LJ diameter, and partial charges, respectively. θ is a bond angle, u eq is the equilibrium value for that bond angle and k u is the force constant.

Groningen molecular simulation (GROMOS) [41]
The force field is defined as: where k b , k u , k f , k j are the force constants, b 0 , u 0 , cos(d), j 0 are the equilibrium values for the respective force constants, q i are the partial charge and C12 ij = 4e ij s 12 ij and C6 ij = 4e ij s 6 ij are the parameters for the Lennard-Jones potential.

Results
Free energy of hydration (ΔG hyd ), free energy of solvation (ΔG solv ) and solubility were determined for 77 molecular systems taken from LigParGen [34][35][36] in order to determine the relationship among these properties. Free energies were calculated through thermodynamic integration while solubility was obtained from density profiles of liquid/liquid simulations of each system mixed with water model TIP4Pe [42]. Solubility is determined by density profiles of binary mixtures with water, in Figure 1 the different cases of solubility are illustrated. These calculations were also performed in the case of ethanol, hexanol and pentane for three different force fields, namely, TraPPE, CHARMM and GROMOS; thus, calculations were performed for 86 ionic liquids finally.
In order to illustrate the different kind of results in Figure 2 the three cases of solubility and their free energy calculations are showed. Figure 2(A-C) shows the density profile of the water/ methanediol mixture, ΔG hyd for methanediol and ΔG solv for methanediol, respectively. Methanediol showed a total solubility with water as can be seen in the density profile, on the other hand, Figure 2(B) presents ΔG hyd with the corresponding λ parameter (red bars) and the accumulative result of this changes (black bars), namely, the integration result; this interpretation is the same in Figure 2(C) for ΔG solv . Figure 2(D-F) shows the results for 2-methylpyridine, which show partial solubility with water and, finally, Figure 2(G-I) presents the results for nitrobenzene, an insoluble molecular system in water. The three molecular systems were taken from LigParGen.
This was the criteria through which results were obtained and analysed for all systems in this work.
As it was mentioned, free energy of hydration (ΔG hyd ) was determined through thermodynamic integration, where the intermolecular interaction parameters of the molecules solvated by water are switched off by perturbing the potential with the λ parameter. In Figure 3 it is showed the result of this property for the mentioned 86 molecular systems. Previously, the solubility in water has been determined and for that reason it was possible to plot ΔG hyd indicating solubility by colours, namely, black for total solubility, red for partial solubility and green corresponds to insolubility. In the figure, a tendency of the property with respect to solubility can not be observed. For this reason, it is not clear how to relate the numerical results of ΔG hyd with solubility. In this figure, also, are included results for three molecules (pentane, ethanol and hexanol) for TraPPE [39,40] (empty circles), GROMOS [41] (empty squares) and CHARMM [37,38] (empty triangles) force field, in order to show that relation among these properties does not depend on the force field.
ΔG solv calculation was also performed for the molecular systems and the numerical results with solubility are shown in   However, by calculating D s =ΔG solv -ΔG hyd as a balance of energies of two systems solvating molecule A (water and molecular system A itself), it is possible to determine a tendency of this result with solubility. In Figure 5 these calculations for the molecular systems mentioned and those from TraPPE, GRO-MOS and CHARMM can be observed. In the figure, it can be seen that as long as the system is soluble in water (black) D s is greater than −6.5 kJ/mol and for insoluble systems (green) D s is less than −6.5 kJ/mol. In the middle, D s for the partially soluble systems fluctuates around −6.5 kJ/mol. For a clear comparison, Figures 3,4 and 5 have the same scale. Tables 1, 2 and 3 in Supporting Information present the organic liquids studied and results of the properties. Table 4 in Supporting Information contains the results for ethanol, hexanol and pentane from TraPPE, CHARMM and GROMOS force fields.

Functional groups
A functional group such as alcohol or ketone normally is by nature, polar, however, as it increases the number of carbon atoms of the hydrocarbon chain attached to it the hydrophobic behaviour increases also its effect. As an example methanol (C1), ethanol (C2) and propanol (C3) are soluble in water while butanol (C4) is partially soluble and hexanol (C6) is basically insoluble. Force fields available are not capable to reproduce this phenomenon since they were not parameterised with the solubility in water or some solution property. However, even though the property is not well reproduced by adding carbon atoms to the functional group the solubility in water decreases just for the interactions. It is not a goal of this work to reproduce the experimental value of some specific property, but also, it is to establish the relation between free energies and solubility.
In Figure 6, it is illustrated the growth of the hydrocarbon chain from dimethylether to butyl-methylether. As it can be seen the hydrocarbon chain grew up until four units. Free energy results of these molecules are included in Figures 3,4  (Colour online) D s defined as ΔG solv -ΔG hyd for LigParGen molecules soluble (black diamonds), partially soluble (red diamonds) and insoluble (green diamonds) in water. Also, circles correspond to three molecules of TraPPE model and squares correspond to GROMOS and triangles to CHARMM force field. The uncertainty in this property is ± 1.19 kJ/mol.  and 5, but, in order to analyse the behaviour of solubility and free energies with the increasing hydrocarbon chain size, in Table 1 the results for ethers, amines, caboxylic acids, thiols and alcohols growing the hydrocarbon chain are shown. As in the case of ethers, all these results have been included in Figures 3,4 and 5. In order to show the tendency of these properties as the hydrophobic structure grows in Figure 7(A-C) it is plotted ΔG hyd , ΔG solv and D s respectively for five different functional groups with an increasing hydrocarbon chain. As an example, amines (black points) are methylamine (C1), ethylamine (C2), propylamine (C3), butylamine (C4), pentylamine (C5) and hexylamine (C6). For alcohols (blue points), carboxylic acids (green points), ethers (red points), and thiols (pink points) the results are also shown in the figures.
For amines, systems are soluble until C5, pentylamine, where ΔG hyd = −18.65 kJ/mol and ΔG solv = −25.02 kJ/mol. For ethers, this transition occurs with methyl-ethylether C2 (ΔG hyd = −8.11 kJ/mol and ΔG solv = −12.32 kJ/mol) and methyl-propylether C3 (ΔG hyd = −6.7 kJ/mol and ΔG solv = −14.86 kJ/mol). On the other hand, for carboxylic acids, the transition takes place with ethanoic acid C2, where ΔG hyd = −40.65 kJ/mol and ΔG solv = −43.22 kJ/mol. In all mentioned cases, it has been observed that in some sizes of the hydrocarbon chain the molecules become insoluble but the free energies do not show a relation among them, even though all mentioned systems describe the same, namely, a functional group which becomes insoluble by adding carbon atoms to the hydrocarbon chain.
However, Figure 7(C) shows D s as the difference between ΔG hyd and ΔG solv . Also, in this figure, it is showed a horizontal line in −6.5 kJ/mol, which is the value plotted in Figure 5 that represents the frontier between soluble and insoluble systems. It is clear that the transition soluble/insoluble for four of the functional groups studied here is around the value, for amines in C5 pentylamine, for alcohols is in C4 butanol and C5 pentanol, for carboxylic acids is C4 butanoic acid and for ethers is in C2 methyl-ethylether and C3 methyl-propylether and for this system, in Figure 7(D), we show the density profiles. Even in the case of thiols, where all systems were insoluble (see Table 1), all D s results are under -6.5 kJ/mol.
It is interesting to observe that in the case of carboxylic acids in ΔG hyd and ΔG solv there is an important change going from C1 to C2, namely, methanoic acid and ethanoic acid respectively (green line in Figure 7(A,B)). In Figure 8 both, C1 and C2 molecules from the functional groups studied here are present, also, atomic electric charges given by the force field OPLS-LigParGen are indicated. Solubility of molecules is a consequence of intermolecular interaction and Coulombic forces play a very important role in this property, therefore the molecular dipole moment was analysed. For carboxylic acids, the difference between the dipole moment of C1 molecule and C2 molecule was found as the biggest among all the functional groups studied, Dm C1−C2 =0.0347 D while in the case of amines Dm C1−C2 =0.0006. This difference indicates that Coulombic interactions are very similar for alcohols and should present important differences for carboxylic acids, as it was observed in ΔG hyd and ΔG solv (Figure 7(A,B)).
Finally, even though the goal of this work is not to reparameterise force fields in order to obtain solubilities according  Figure 5). Colours correspond to: black/amines, red/ethers, green/carboxylic acids, blue/alcohols and pink/thiols. (d) Density profiles of the ethers/water mixtures, colours correspond to: blue to dimethylether, black to ethyl-methylether, red to propyl-methylether and green to butyl-methylether.
to experimental data, it is interesting to show how simulation results correlate with experiments, as long as the experimental result was available. In Figure 9 the correlation between experimental data simulation results can be observed. Numerical data and molecules in this figure are described in the Supporting Information.
This figure shows that solubility in water, in general, is not well reproduced by the force fields used in this work. This result has been exposed before [15,29] and represents an area of opportunity for the improvement of the parameterisation methods

Conclusions and discussion
Eighty-six organic liquids were studied through free energy calculations and solubility in water. ΔG hyd does not present a correlation with solubility in water and that is the case also for ΔG solv . In a liquid/liquid system formed of water and an organic liquid A, a single molecule A in the interface, since behaviour is a dynamical process, has two options, being solvated by water or being solvated by system A itself, the final result depends on the interactions of such molecule with water and with its molecular system A. For that reason it is not enough to evaluate ΔG hyd to find out if system A is soluble in water, instead, it has to be also determined ΔG solv and make a comparison between these free energies. That is how D s is defined as ΔG solv -ΔG hyd .
It was found a numerical relation between D s and solubility, namely, insoluble systems present a D s ≪6.5 kJ/mol and soluble systems a D s ≫ −6.5 kJ/mol, D s for systems which present a partial solubility is around −6.5 kJ/mol.
In order to show that results do not depend on the force field chosen, these simulations were also performed for ethanol, hexanol and pentane of three different force fields (TraPPE, CHARMM and GROMOS), and again, ΔG hyd and ΔG solv did not show a correlation with solubility in water and D s did.
In a molecular structure, as long as the hydrocarbon chain attached to a functional group grows, the hydrophobicity effect increases, for example, in the case of amines, methylamine, ethylamine, propylamine, butylamine, and pentylamine are soluble in water while hexylamine is practically insoluble [43][44][45][46]. This phenomenon occurs in each functional group.
In the case of molecular simulations, independently of the force field, it is a fact that as the hydrocarbon chain grows in a functional group, the hydrophobicity also increases, but the experimental value of solubility, normally, is not well reproduced by force fields available. It is not a goal of this work to reproduce the experimental results, but also, to find a relation between free energy and solubility in water since free energy calculations through MD are fast to obtain in comparison with liquid/liquid simulations.
Liquid/liquid simulations with water and free energy calculations were carried out for amines, ethers, carboxylic acids, alcohols, and thiols to observe the behaviour of these properties as long as the hydrocarbon chain grew. Results from the  simulations verified the non-correlated values between ΔG solv and ΔG hyd with solubility in water and, on the other hand, the mentioned correlation between D s and solubility was obtained again.
This result is valuable because it describes a straight relation between free energy and solubility and allows us to determine, qualitatively, the solubility by free energy calculations (thermodynamic integration of a relatively minor number of molecules), which are much cheaper than those of solubility (semi-isotropic simulations with a considerably major number of molecules).

Disclosure statement
No potential conflict of interest was reported by the author(s).