Coupled cluster theory with the polarizable continuum model of solvation

The environment may significantly affect molecular properties. Thus, it is desirable to account explicitly for these effects on the wave function and its derivatives, especially when the latter are evaluated with accurate methods, such as those belonging to coupled cluster (CC) theory. In this tutorial review, we discuss how to combine CC methods with the polarizable continuum model of solvation (PCM). We describe useful approximations that include the solvent response to the correlation and excited state equations while maintaining the computational cost comparable to in vacuo calculations. Although applied to PCM, the theoretical framework presented in this review is general and can be used with any polarizable embedding model. Representative applications of the CC-PCM method to ground and excited state properties of solvated molecules are presented, and comparisons with experiment, and between the full and approximate schemes are discussed.     

Molecular orientations at interfaces by extended polarizable continuum model

This study presents an investigation into orientation of molecular solutes at the interface of liquid water and other media. The calculation of electrostatic free energy of molecular solute is based on an extension of the polarizable continuum model (PCM) to interfacial system. The extended PCM computational scheme is incorporated with the self‐consistent field procedure which is necessary to obtain more accurate electrostatic free energy and charge density distribution. The computation of non‐electrostatic energy for interfacial system is also realized. Applying the numerical procedure to molecular systems, N,N′‐diethyl‐p‐nitroaniline (DEPNA) at air/water interface and p‐nitrophenol (PNP) at cyclohexane/water interface, the average orientational angles are in reasonable agreement with the experimental results. Taking both the electrostatic and the non‐electrostatic energies into account, the analysis on the energy profiles shows that the electrostatic solvation energy is the dominant factor in determining the orientation angle for PNP, whereas for DEPNA, the orientation angle mainly depends on the cavitation energy. This suggests that, in addition to the electrostatic energy, taking the cavitation energy into account may provide a more complete view when we survey the molecular orientation at interface. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Achieving linear-scaling computational cost for the polarizable continuum model of solvation

This work describes a new and low-scaling implementation of the polarizable continuum model (PCM) for computing the self-consistent solvent reaction field. The PCM approach is both general and accurate. It is applicable in the framework of both quantum and classical calculations, and also to hybrid quantum/classical methods. In order to further extend the range of applicability of PCM we addressed the problem of its computational cost. The generation of the finite-elements molecular cavity has been reviewed and reimplemented, achieving linear scaling for systems containing up to 500 atoms. Linear scaling behavior has been achieved also for the iterative solution of the PCM equations, by exploiting the fast multipole method (FMM) for computing electrostatic interactions. Numerical results for large (both linear and globular) chemical systems are discussed.     

Parallelization of the integral equation formulation of the polarizable continuum model for higher-order response functions

We present a parallel implementation of the integral equation formalism of the polarizable continuum model for Hartree-Fock and density functional theory calculations of energies and linear, quadratic, and cubic response functions. The contributions to the free energy of the solute due to the polarizable continuum have been implemented using a master-slave approach with load balancing to ensure good scalability also on parallel machines with a slow interconnect. We demonstrate the good scaling behavior of the code through calculations of Hartree-Fock energies and linear, quadratic, and cubic response function for a modest-sized sample molecule. We also explore the behavior of the parallelization of the integral equation formulation of the polarizable continuum model code when used in conjunction with a recent scheme for the storage of two-electron integrals in the memory of the different slaves in order to achieve superlinear scaling in the parallel calculations.     

Alkaline hydrolysis of ethylene phosphate: an ab initio study by supermolecule model and polarizable continuum approach

The alkaline hydrolysis reaction of ethylene phosphate (EP) has been investigated using a supermolecule model, in which several explicit water molecules are included. The structures and single-point energies for all of the stationary points are calculated in the gas phase and in solution at the B3LYP/6-31++G(df,p) and MP2/6-311++G(df,2p) levels. The effect of water bulk solvent is introduced by the polarizable continuum model (PCM). Water attack and hydroxide attack pathways are taken into account for the alkaline hydrolysis of EP. An associative mechanism is observed for both of the two pathways with a kinetically insignificant intermediate. The water attack pathway involves a water molecule attacking and a proton transfer from the attacking water to the hydroxide in the first step, followed by an endocyclic bond cleavage to the leaving group. While in the first step of the hydroxide attack pathway the nucleophile is the hydroxide anion. The calculated barriers in aqueous solution for the water attack and hydroxide attack pathways are all about 22 kcal/mol. The excellent agreement between the calculated and observed values demonstrates that both of the two pathways are possible for the alkaline hydrolysis of EP.     

Solvation effects on reaction profiles by the polarizable continuum model coupled with the Gaussian density functional method

An efficient version of the polarizable continuum model for solvation has been implemented in the Gaussian density-functional-based code called deMon. Solvation free energies of representative compounds have been calculated as a preliminary test. The hydration effects on the reaction profile of the Cl⁻+CH₃Cl→ClCH₃+Cl⁻ reaction and the thermodynamics of the Menschutkin reaction have also been investigated. Finally, the conformational behavior of the 1,2-diazene cis–trans isomerization process in water was examined. Comparisons between the results obtained and the available experimental data and previous theoretical computations have been made. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 290–299, 1998     

Polarizable continuum model with the fragment molecular orbital-based time-dependent density functional theory

The polarizable continuum model (PCM) for describing the solvent effect was combined with the fragment molecular orbital-based time-dependent density functional theory (TDDFT). Several levels of the many-body expansion were implemented, and the importance of the many-body contributions to the singlet-excited states was discussed. To calibrate the accuracy, we performed a number of the model calculations using our method and the regular TDDFT in solution, applying them to phenol and polypeptides at the long-range corrected BLYP/6-31G* level. It was found that for systems up to 192 atoms the largest error in the excitation energy was 0.006 eV (vs. the regular TDDFT/PCM of the full system). The solvent shifts and the conformer effects were discussed, and the scaling was found to be nearly linear. Finally, we applied our method to the lowest singlet excitation of the photoactive yellow protein (PYP) in aqueous solution and determined the excitation energy to be in reasonable agreement with experiment. The excitation energy analysis provided the contributions of individual residues, and the main factors as well as their solvent shifts were determined.     

An improved iterative solution to solve the electrostatic problem in the polarizable continuum model

 We present a discrete iterative interpolation scheme (DIIS) to improve the convergence rate of electrostatic calculations in the polarizable continuum model (PCM) to describe solvent effects on molecular solutes. The electrostatic calculations may easily become the bottleneck of the calculation when the solute size is large. For large molecules iterative procedures turn out to be computationally more convenient than matrix inversion or closure methods. The DIIS scheme is compared here to another iterative procedure (DAMP) and to the biconjugate gradient (BCG) method. The comparisons show that DIIS leads to a sizeable saving of computational time for the C-PCM and IEF-PCM methods (average 40%) compared to DAMP, and more than 50% with respect to the BCG method. Received: 5 October 2000 / Accepted: 13 November 2000 / Published online: 19 January 2001     

Implicit solvent models: Combining an analytical formulation of continuum electrostatics with simple models of the hydrophobic effect

The recent development of approximate analytical formulations of continuum electrostatics opens the possibility of efficient and accurate implicit solvent models for biomolecular simulations. One such formulation (ACE, Schaefer & Karplus, J. Phys. Chem., 1996, 100:1578) is used to compute the electrostatic contribution to solvation and conformational free energies of a set of small solutes and three proteins. Results are compared to finite-difference solutions of the Poisson equation (FDPB) and explicit solvent simulations and experimental data where available. Small molecule solvation free energies agree with FDPB within 1–1.5 kcal/mol, which is comparable to differences in FDPB due to different surface treatments or different force field parameterizations. Side chain conformation free energies of aspartate and asparagine are in qualitative agreement with explicit solvent simulations, while 74 conformations of a surface loop in the protein Ras are accurately ranked compared to FDPB. Preliminary results for solvation free energies of small alkane and polar solutes suggest that a recent Gaussian model could be used in combination with analytical continuum electrostatics to treat nonpolar interactions. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 322–335, 1999     

Using polarizable POSSIM force field and fuzzy‐border continuum solvent model to calculate pKa shifts of protein residues

Our Fuzzy‐Border (FB) continuum solvent model has been extended and modified to produce hydration parameters for small molecules using POlarizable Simulations Second‐order Interaction Model (POSSIM) framework with an average error of 0.136 kcal/mol. It was then used to compute pK _a shifts for carboxylic and basic residues of the turkey ovomucoid third domain (OMTKY3) protein. The average unsigned errors in the acid and base pK _a values were 0.37 and 0.4 pH units, respectively, versus 0.58 and 0.7 pH units as calculated with a previous version of polarizable protein force field and Poisson Boltzmann continuum solvent. This POSSIM/FB result is produced with explicit refitting of the hydration parameters to the pK _a values of the carboxylic and basic residues of the OMTKY3 protein; thus, the values of the acidity constants can be viewed as additional fitting target data. In addition to calculating pK _a shifts for the OMTKY3 residues, we have studied aspartic acid residues of Rnase Sa. This was done without any further refitting of the parameters and agreement with the experimental pK _a values is within an average unsigned error of 0.65 pH units. This result included the Asp79 residue that is buried and thus has a high experimental pK _a value of 7.37 units. Thus, the presented model is capable or reproducing pK _a results for residues in an environment that is significantly different from the solvated protein surface used in the fitting. Therefore, the POSSIM force field and the FB continuum solvent parameters have been demonstrated to be sufficiently robust and transferable. © 2016 Wiley Periodicals, Inc.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

Geometry optimization of molecular structures in solution by the polarizable continuum model

A new implementation of analytical gradients for the polarizable continuum model is presented, which allows Hartree-Fock and density functional calculations taking into account both electrostatic and nonelectrostatic contributions to energies and gradients for closed and open shell systems. Simplified procedures neglecting the derivatives of the cavity surface and/or using single spheres for XH_n groups have also been implemented and tested. The solvent-induced geometry relaxation has been studied for a number of representative systems in order to test the efficiency of the procedure and to investigate the role of different contributions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 404–417, 1998     

Coupled‐cluster theory for the polarizable continuum model. III. A response theory for molecules in solution

A coupled‐cluster (CC) response functions theory for molecular solutes described with the framework of the polarizable continuum model (PCM) is presented. The theory is an extension to the dynamical molecular properties of the PCM‐CC analytic derivatives recently proposed for the calculation of static molecular properties (Cammi, Jr Chem Phys 2009, 131, 164104). The theory is presented for linear and quadratic response functions, and the operative expressions of these response functions can accurately account for the nonequilibrium solvation effects. The excitation energies and transition moments of the solvated chromophores have been determined from the linear response functions. Accurate expressions for gradients of excitation energies for the evaluation of the excited state properties have been also discussed. © 2012 Wiley Periodicals, Inc. Int J Quantum Chem, 2012     

A study of amino-protecting groups using the polarizable continuum model (PCM)

In order to better understand the performance of 1,2-dimethyl-5-acetyl barbituric acid (DMB) as an amino protecting group relative to 5,5-dimethylcyclohexane-1,3-dione (DMD), ab initio calculations were performed. pK_a calculations using the PCM model indicated that both molecules are more acidic in the enol form. Therefore, the protecting reaction of these molecules should involve the anions formed from the loss of a proton from the enol compounds. Contrary to what would be expected, the larger efficiency exhibited by the DMB molecule cannot be attributed to an extension of the electronic conjugation effect. In the absence of any other noticeable effect that could be responsible for the greater efficiency of the DMB molecule, we are inclined to believe that the difference could be accounted for by the presence of two independent centers of conjugation.This paper is dedicated to Jacopo Tomasi in recognition of his outstanding contribution to the field of computational chemistry in solution. The authors are honored to contribute to this volume; especially so for two of them (COS and MACN) who have the privilege of his friendship.Acknowledgements The authors would like to thank the Brazilian research agencies CNPq, CAPES and FAPERJ for the financial support. C. O. da Silva thanks the Dipartimento di Chimica e Chimica Industriale, University of Pisa, where the MCSCF calculations were performed.     

A study on orientation and absorption spectrum of interfacial molecules by using continuum model

In this work, a numerical procedure based on the continuum model is developed and applied to the solvation energy for ground state and the spectral shift against the position and the orientation of the interfacial molecule. The interface is described as a sharp boundary separating two bulk media. The polarizable continuum model (PCM) allows us to account for both electrostatic and nonelectrostatic solute-solvent interactions when we calculate the solvation energy. In this work we extend PCM to the interfacial system and the information about the position and orientation of the interfacial molecule can be obtained. Based on the developed expression of the electrostatic free energy of a nonequilibrium state, the numerical procedure has been implemented and used to deal with a series of test molecules. The time-dependent density functional theory (TDDFT) associated with PCM is used for the electron structure and the spectroscopy calculations of the test molecules in homogeneous solvents. With the charge distribution of the ground and excited states, the position- and orientation-dependencies of the solvation energy and the spectrum have been investigated for the interfacial systems, taking the electrostatic interaction, the cavitation energy, and the dispersion-repulsion interaction into account. The cavitation energy is paid particular attention, since the interface portion cut off by the occupation of the interfacial molecule contributes an extra part to the stabilization for the interfacial system. The embedding depth, the favorable orientational angle, and the spectral shift for the interfacial molecule have been investigated in detail. From the solvation energy calculations, an explanation has been given on why the interfacial molecule, even if symmetrical in structure, tends to take a tilting manner, rather than perpendicular to the interface.