Molecular dynamics simulations have been performed to compute the solvation free energy and the octanol/water partition coefficients for a challenging set of selected organic molecules, characterized by the simultaneous presence of functional groups coarsely spanning a large portion of the chemical space in drug-like compounds and, in many cases, by a complex conformational landscape (2-propoxyethanol, acetylsalicylic acid, cyclohexanamine, dialifor, ketoprofen, nitralin, profluralin, terbacil). OPLS-AA and GAFF2 parameterizations of the organic molecules and of 1-octanol have been done via the web-based automatic parameter generators, LigParGen [Dodda et al.

In this work, we investigate the effects of different atomic charge and nonelectrostatic models on the hydration energies of neutral molecules, using an implicit solvation model. The solvation free energy is divided into two main components, the first resulting from a self‐consistent reaction field treatment of the bulk electrostatics obtained by solving the Poisson equation in a finite‐difference (FD) approach where the solute charge density is approximated by atomic charges, the second corresponding to short‐range interactions between the solute and the solvent in the first solvation shell. Five different atomic charge models (Mulliken, Hirshfeld, Hirshfeld‐I, CM5 and its iterative version, CM5‐I) have been considered, both at the Hartree–Fock (HF) and B3LYP levels, with three different basis sets, alongside two nonelectrostatic models including the cavity, dispersion, and solvent structural effects (CDS) model. Averaging over the three considered basis sets, Hirshfeld charges combined to the CDS model led to the lowest mean unsigned error (MUE), with a value of 0.92 kcal/mol with respect to the experimental data. On the other hand, a MUE of 2.02 kcal/mol was obtained with CM5 charges combined to the CDS model, highlighting the low transferability of the original CDS parameters developed for the generalized Born electrostatics to a different electrostatics model. By scaling down the CM5 charges to better balance with the original CDS model, a MUE of 0.68 kcal/mol was however obtained, outlining the delicate balance existing between the electrostatic and nonelectrostatic contributions to the solvation free energy in implicit solvation models.

The structural and dynamic properties of imidazole in aqueous solution have been studied by means of classical and ab initio molecular dynamics simulations. We developed a new force field for the imidazole molecule with improved modelling of the electrostatic interactions, specifically tailored to address the well known drawbacks of existing force fields based on the atomic fractional charges approach. To this end, we reparametrized the charge distribution on the heterocyclic ring, introducing an extra site accounting for the lone pair on the deprotonated nitrogen. The accuracy of the model in describing the hydrogen bond pattern in aqueous solvent has been confirmed by comparing the classical results on imidazole-water interactions to accurate Car-Parrinello molecular dynamics simulations. The proposed classical model reproduces satisfactorily the experimental water/octanol partition coefficient of imidazole, as well as the structure of the imidazole molecular crystal. The force field has been finally applied to simulate aqueous solutions at various imidazole concentrations to obtain information on both imidazole-water and imidazole-imidazole interactions, providing a description of the different molecular arrangements in solution.

We present a comparison of the performances of three continuum solvation models, namely the Solvation Model Density (SMD), VASPsol and Finite-Difference Poisson Boltzmann (FDPB) approaches, for the calculation of the hydration energies of molecules, polymers and semiconductor surfaces. For finite molecular systems, all three models have been considered, and the computed data have been compared to available experimental solvation energies for a test set of 630 neutral solutes. For infinite periodic systems, due to the lack of periodic implementations of the SMD model and of experimental solvation energies, only a comparison between the data obtained with the VASPsol and FDPB approaches has been performed. As a key criterion of the validation of a periodic implementation, the size-extensivity of the solvation energy of a model of a poly glycine chain has been considered. In addition, the effect of the surface orientation and of the slab thickness on the computed solvation energy has been investigated by considering five low-index surface orientations of galena PbS, due to the importance of this material in environmental-related processes. For finite molecular systems, all models performed satisfactorily on the whole test set with mean unsigned errors (MUE) values close to the target 1 kcal/mol error threshold with respect to the experimental data. However, depending on the chosen solvation model, some solute classes revealed challenging with MUE reaching 4.5 kcal/mol in the worst cases. Overall, the best performing models were found to be: FDPB > SMD ≈ VASPsol. For infinite periodic systems, the size-extensivity of the computed solvation energies has been verified, both with the VASPsol and FDPB models. In addition, calculations performed on PbS surface models further revealed at least a qualitative agreement between the two solvation models. In particular, the trends in surface energy variation between gasphase and solvated cases, computed solvation energies, as well as the surface solvation energy convergence as a function of the number of layers in the slab models have all been found very similar between these two approaches. Overall, therefore, these results are encouraging for the applicability of continuum solvation models to a broad range of solutes, ranging from finite small molecular solutes to extended periodic systems.

In this article, we explore an alternative to the analytical Gauss-Bonnet approach for computing the solvent-accessible surface area (SASA) and its nuclear gradients. These two key quantities are required to evaluate the nonelectrostatic contribution to the solvation energy and its nuclear gradients in implicit solvation models. We extend a previously proposed analytical approach for finite systems based on the stereographic projection technique to infinite periodic systems such as polymers, nanotubes, helices, or surfaces and detail its implementation in the Crystal code. We provide the full derivation of the SASA nuclear gradients, and introduce an iterative perturbation scheme of the atomic coordinates to stabilize the gradients calculation for certain difficult symmetric systems. An excellent agreement of computed SASA with reference analytical values is found for finite systems, while the SASA size-extensivity is verified for infinite periodic systems. In addition, correctness of the analytical gradients is confirmed by the excellent agreement obtained with numerical gradients and by the translational invariance achieved, both for finite and infinite periodic systems. Overall therefore, the stereographic projection approach appears as a general, simple, and efficient technique to compute the key quantities required for the calculation of the nonelectrostatic contribution to the solvation energy and its nuclear gradients in implicit solvation models applicable to both finite and infinite periodic systems. K E Y W O R D S crystal, implicit solvation, solvent-accessible surface area