Essential Computational Modeling in Chemistry

By Philippe G. Ciarlet

Essential Computational Modeling in Chemistry presents key contributions selected from the volume in the Handbook of Numerical Analysis: Computational Modeling in Chemistry Vol. 10(2005).

Computational Modeling is an active field of scientific computing at the crossroads between Physics, Chemistry, Applied Mathematics and Computer Science. Sophisticated mathematical models are increasingly complex and extensive computer simulations are on the rise. Numerical Analysis and scientific software have emerged as essential steps for validating mathematical models and simulations based on these models. This guide provides a quick reference of computational methods for use in understanding chemical reactions and how to control them. By demonstrating various computational methods in research, scientists can predict such things as molecular properties. The reference offers a number of techniques and the numerical analysis needed to perform rigorously founded computations.

Various viewpoints of methods and applications are available for researchers to chose and experiment with; Numerical analysis and open problems is useful for experimentation; Most commonly used models and techniques for the molecular case is quickly accessible

• Language: English
• Publisher: Elsevier Science
• Release date: Dec 7, 2010
• ISBN: 9780444537614

Book preview

Essential Computational Modeling in Chemistry

A derivative of Handbook of Numerical Analysis Special Volume: Computational Chemistry, Vol 10

First Edition

P.G. Ciarlet

Laboratoire Jacques-Louis Lions, Université Pierre et Marie Curie, 4 Place Jussieu, 75005 PARIS, France

Department of Mathematics, City University of Hong Kong, Tat Chee Avenue, KOWLOON, Hong Kong

C. Le Bris

CERMICS, Ecole Nationale des Ponts et Chaussées, 77455 Marne La Vallée, France

ELSEVIER

Amsterdam • Boston • Heidelberg • London • New York • Oxford • Paris

San Diego • San Francisco • Singapore • Sydney • Tokyo

Table of Contents

Cover image

Title page

Copyright page

General Preface

Contributors

The Modeling and Simulation of the Liquid Phase

Chapter I Introduction to Liquids

Chapter II Continuum Models

Chapter III Computer Simulations

Chapter IV Hybrid Methods

Chapter V Computing Properties

Chapter VI Excited States

Computational Approaches of Relativistic Models in Quantum Chemistry

1 Introduction

2 Linear Dirac equations

3 The MCDF method for atoms

4 Numerical relativistic methods for molecules

Quantum Monte Carlo Methods for the Solution of the Schrödinger Equation for Molecular Systems

Preface

Part I Introduction

Part II Algorithms

Part III Special topics

Acknowledgments

Finite Difference Methods for Ab Initio Electronic Structure and Quantum Transport Calculations of Nanostructures

1 Introduction

2 Electronic structure calculation by finite differences

3 Quantum transport

4 Applications: conductivity from ab initio local orbital Hamiltonian

Acknowledgments

Simulating Chemical Reactions in Complex Systems

1 Introduction

2 Classical theories of reaction rates

3 Calculating condensed-phase potential energy surfaces

4 Simulation methods for investigating chemical reactions

5 Quantum algorithms

6 Challenges and perspectives

Acknowledgments

Biomolecular Conformations Can Be Identified as Metastable Sets of Molecular Dynamics

1 Introduction

2 Conceptual preliminaries

3 Description of dynamical behavior

4 Metastability

5 Transfer operators

6 Numerical realization

7 Illustrative numerical experiments

8 Application to biomolecular systems

9 Appendix

Numerical Methods for Molecular Time-Dependent Schrödinger Equations – Bridging the Perturbative to Nonperturbative Regime

1 Introduction

2 Gauges and representations

3 Numerical schemes

4 Boundary conditions and energy spectra

5 Beyond the dipole approximation

6 Conclusion

Acknowledgment

Control of Quantum Dynamics: Concepts, Procedures, and Future Prospects

1 Introduction

2 Basic principles

3 Controllability of quantum mechanical systems

4 Quantum control algorithms

5 Challenges for the future

6 Conclusion

Acknowledgments

Subject Index

Copyright

Elsevier

30 Corporate Drive, Suite 400, Burlington, MA 01803, USA

Linacre House, Jordan Hill, Oxford 0X2 8DP, UK

Radarweg 29, PO Box 211, 1000 AE Amsterdam, The Netherlands

First edition 2011

© 2011 Elsevier B.V. All rights reserved

Material in the text originally appeared in the Handbook of Numerical Analysis Special Volume: Computational Chemistry, Vol 10 edited by P.G. Ciarlet (editor) and C. Le Bris (guest editor) (Elsevier B.V., 2003)

No part of this publication may be reproduced, stored in a retrieval system or transmitted in any form or by any means electronic, mechanical, photocopying, recording or otherwise without the prior written permission of the publisher Permissions may be sought directly from Elsevier’s Science & Technology Rights Department in Oxford, UK: phone (+44) (0) 1865 843830; fax (+44) (0) 1865 853333; email: [e-mail]. Alternatively you can submit your request online by visiting the Elsevier web site at http://elsevier.com/locate/permissions, and selecting Obtaining permission to use Elsevier material.

Notice

No responsibility is assumed by the publisher for any injury and/or damage to persons or property as a matter of products liability, negligence or otherwise, or from any use or operation of any methods, products, instructions or ideas contained in the material herein. Because of rapid advances in the medical sciences, in particular, independent verification of diagnoses and drug dosages should be made.

British Library Cataloguing in Publication Data

A catalogue record for this book is available from the British Library

Library of Congress Cataloging-in-Publication Data

A catalog record for this book is available from the Library of Congress

For information on all North Holland publications visit our web site at books.elsevier.com

Printed and bound in Great Britain

11 12 13 10 9 8 7 6 5 4 3 2 1

ISBN: 978–0-444–53754-6

General Preface

In the early eighties, when Jacques-Louis Lions and I considered the idea of a Handbook of Numerical Analysis, we carefully laid out specific objectives, outlined in the following excerpts from the General Preface, which has appeared at the beginning of each of the volumes published so far:

During the past decades, giant needs for ever more sophisticated mathematical models and increasingly complex and extensive computer simulations have arisen. In this fashion, two indissociable activities, mathematical modeling and computer simulation, have gained a major status in all aspects of science, technology and industry.

In order that these two sciences be established on the safest possible grounds, mathematical rigor is indispensable. For this reason, two companion sciences, Numerical Analysis and Scientific Software, have emerged as essential steps for validating the mathematical models and the computer simulations that are based on them.

Numerical Analysis is here understood as the part of Mathematics that describes and analyzes all the numerical schemes that are used on computers; its objective consists in obtaining a clear, precise and faithful representation of all the information contained in a mathematical model; as such, it is the natural extension of more classical tools, such as analytic solutions, special transforms, functional analysis, as well as stability and asymptotic analysis.

The various volumes comprising the Handbook of Numerical Analysis thoroughly cover all the major aspects of Numerical Analysis, by presenting accessible and in-depth surveys, which include the most recent trends.

More precisely, the Handbook covers the basic methods of Numerical Analysis, gathered under the following general headings:

n,

– Finite Difference Methods,

– Finite Element Methods,

– Techniques of Scientific Computing.

It also covers the numerical solution of actual problems of contemporary interest in Applied Mathematics, gathered under the following general headings:

– Numerical Methods for Fluids,

– Numerical Methods for Solids.

In retrospect, it can be safely asserted that Volumes I to IX, which were edited by both of us, fulfilled most of these objectives, thanks to the eminence of the authors and the quality of their contributions.

After Jacques-Louis Lions’ tragic loss in 2001, it became clear that Volume IX would be the last one of the type published hitherto, that is, edited by both of us and devoted to some of the general headings defined above. It was then decided, in consultation with the publisher, that each future volume will instead be devoted to a single "specific application" and called for this reason a "Special Volume". "Specific applications includes Mathematical Finance, Meteorology, Celestial Mechanics, Computational Chemistry, Living Systems, Electromagnetism, Computational Mathematics, etc. It is worth noting that the inclusion of such specific applications" in the Handbook of Numerical Analysis was part of our initial project.

To ensure the continuity of this enterprise, I will continue to act as Editor of each Special Volume, whose conception will be jointly coordinated and supervised by a Guest Editor.

P.G. Clarlet, July 2002

Contributors

Alán Aspuru-Guzik; William A. Lester, Jr.     Department of Chemistry, University of California, Berkeley, California, USA

André D. Bandrauk; Hui-Zhong Lu     Laboratoire de Chimie Théorique, Faculté des Sciences, Université de Sherbrooke, Canada

Eric Brown, [e-mail]     Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ, USA. (E. Brown)

J.P. Desclaux, [e-mail]     Sassenage, France.

J. Dolbeault, [e-mail]     CEREMADE, Unité Mixte de Recherche du CNRS no. 7534 et Université Paris IX- Dauphine, Paris cedex, France.

M.J. Esteban, [e-mail]     CEREMADE, Unité Mixte de Recherche du CNRS no. 7534 et Université Paris IX- Dauphine, Paris cedex, France.

Jean-Luc Fattebert, [e-mail]     Center for Applied Scientific Computing (CASC), Lawrence Livermore National Laboratory, P.O. Box 808, L-561, Livermore, CA, 94551, USA

M.J. Field     Laboratoire de Dynamique Moléculaire, Institut de Biologie Structurale – Jean-Pierre Ebel, 41 rue Jules Horowitz, Grenoble cedex 1, France

Wilhelm Huisinga, [e-mail]     Institute of Mathematics II, Department of Mathematics and Computer Science, Free University (FU) Berlin, Germany

P. Indelicato, [e-mail]     Laboratoire Kastler-Brossel, Unité Mixte de Recherche du CNRS no. C8552, École Normale Supérieure et Université Pierre et Marie Curie, Paris cedex, France.

Patrick Laug     Dipartimento di Chimica e Chim. Indus., Via Risorgimento 35, University of Pisa, Pisa, Italy INRIA Rocquencourt, Gamma project, BP 105, Le Chesnay cedex, France

B. Mennucci     Dipartimento di Chimica e Chim. Indus., Via Risorgimento 35, University of Pisa, Pisa, Italy INRIA Rocquencourt, Gamma project, BP 105, Le Chesnay cedex, France

Marco Buongiorno Nardelli     Department of Physics, North Carolina State University, Raleigh, NC, and Center for Computational Sciences (CCS) and Computational Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN, 37830, USA

Herschel Rabitz, [e-mail]     Department of Chemistry, Princeton University, Princeton, NJ, USA.

Christof Schütte, [e-mail]     Institute of Mathematics II, Department of Mathematics and Computer Science, Free University (FU) Berlin, Germany

E. Séré, [e-mail]     CEREMADE, Unité Mixte de Recherche du CNRS no. 7534 et Université Paris IX- Dauphine, Paris cedex, France.

J. Tornasi     Dipartimento di Chimica e Chim. Indus., Via Risorgimento 35, University of Pisa, Pisa, Italy INRIA Rocquencourt, Gamma project, BP 105, Le Chesnay cedex, France

Gabriel Turinici, [e-mail], [e-mail]     INRIA Rocquencourt, B.P. 105, Le Chesnay cedex, France CERMICS-ENPC, Champs sur Marne, Marne la Vallée cedex, France. (G. Turinici)

The Modeling and Simulation of the Liquid Phase

J. Tomasi    Dipartimento di Chimica e Chim. Indus., Via Risorgimento 35, University of Pisa, Pisa, Italy

INRIA Rocquencourt, Gamma project, BP 105, Le Chesnay cedex, France

B. Mennucci    Dipartimento di Chimica e Chim. Indus., Via Risorgimento 35, University of Pisa, Pisa, Italy

INRIA Rocquencourt, Gamma project, BP 105, Le Chesnay cedex, France

Patrick Laug    Dipartimento di Chimica e Chim. Indus., Via Risorgimento 35, University of Pisa, Pisa, Italy

INRIA Rocquencourt, Gamma project, BP 105, Le Chesnay cedex, France

Contents

CHAPTER I INTRODUCTION TO LIQUIDS   5

1. Physical approaches to the study of liquids   6

2. Chemical approaches to the study of liquids   11

CHAPTER II CONTINUUM MODELS   16

3. The energy of the system and the solvation energy   20

4. Computation of the electrostatic contribution   26

5. Geometry optimization   33

6. Molecular surface meshing   34

CHAPTER III COMPUTER SIMULATIONS   52

7. Monte Carlo methods   52

8. Molecular Dynamics   55

CHAPTER IV HYBRID METHODS   58

9. Quantum mechanics/molecular mechanics (QM/MM)   58

10. Layered methods: ONIOM   60

11. The effective fragment potential (EFP) method   61

12. Car-Parrinello ab initio molecular dynamics: AIMD   62

CHAPTER V COMPUTING PROPERTIES   65

13. Liquids   65

14. Molecules in solution   71

CHAPTER VI EXCITED STATES   90

15. Continuum solvation approaches   91

REFERENCES   95

An overview on the theoretical and computational methodologies developed so far to study liquids and solutions is presented. The main characteristics of the different methods are outlined and the advantages and shortcomings of the computational approaches are discussed. Particular attention is focused on a specific class of methods, known as continuum solvation models, for which a more detailed description of theoretical and computational aspects is presented. For this class of methods, the concept of molecular cavity is introduced together with a description of the numerical techniques developed to mesh the corresponding surface. For selected methods representing those of larger use, an overview on their applications to the evaluation of energies and properties of liquid systems and molecules in solution is also presented.

Chapter I Introduction to Liquids

For very long tradition chemistry has to deal with liquids. The chain of chemical manipulations performed in laboratories and in factories are mostly performed in a liquid medium and one of the first questions asked when a new chemical is introduced in the use is what is its solubility in solvents. There is a large number of solvents in use in the chemical practice, more than two thousand, and this number is not sufficient: often it is convenient to use solvent mixtures. Chemical reactions, the core of chemistry, are very sensitive to the solvent used to put in intimate contact the reagent species: often a delicate reaction fails if the solvent was not properly chosen. Important areas of chemical research are inherently tied to a solvent. An important example is that of biological systems. All the complex machinery of biological systems works only if the system is in water: without this solvent event the intimate structures of biological molecules collapse and loose their activity.

Chemists are so accustomed to consider many properties of liquids, from the basic physical properties to others, more specific and playing a role in specific cases in the large variety of problems that chemistry has to face. We cannot give here an overview of the properties of interest for liquids, many among which could be the object of computer modeling and we limit ourselves to indicate a few points more directly related to what will be exposed in the following of this chapter.

Liquids often occur in the chemical practice in large quantities. The effects due to the spatial limitations of the liquid sample have not great influence on the phenomena occurring at the interior, which in the current practice is called the bulk solvent. On large scale the bulk pure liquids and solutions exhibit isotropy. Liquids can be however dispersed. A typical example are the water droplets dispersed in the atmosphere, as for example in the fog. Other examples occur in liquid mixtures that present a miscibility lacuna. In all cases the ratio surface region/bulk can be relatively large and phenomena occurring at the interface may present properties different with respect to the same occurring in the bulk.

The region of separation between a bulk liquid and a solid body has in some cases a decisive importance. One example is the separation between a salt solution and a metal electrode. All the electrochemical phenomena are influenced by the behavior of this thin portion of the liquid. Another example are the phenomena leading to dissolution of solid bodies in the liquid phase, and to the deposition of material components of the liquid into solid particles. The surface of separation of a liquid and a solid surface presents features of different type when the liquid is in a limited amount. The liquid phase can be present under the form of separate drops, or of a thin liquid film on the solid surface: both occurrences have a great importance in some chemical problems.

Surfaces of bulk liquids can be covered by a thin layer of an immiscible liquid. This is a phenomenon frequently observed in everyday life. Such thin layers can be organized into a structure that maintains some aspects of normal liquids but presents a high local order. With an appropriate selection of the chemical composition of this second phase, it is easy to form ordered layers of well-defined molecular thickness and membranes. Membranes, on their part, preserves fluid-like properties at their interior which can be exploited, for example, in the machinery of biological systems. A liquid can contain small solid particles within the bulk. Many phenomena and many practical applications are based on this particular type of liquid systems. A well-known phenomenon that has played an important role in the development of science is the Brownian motion. The solvent eases the dispersion in the bulk of such small particles, that are however subjected to gravity forces leading to sedimentation. Gravity, one of the basic forces in the universe, plays a little role in the microscopic approach to the study of material systems, but this example shows that there are reasons to forecast that in the future also gravity will be included into some computational models. A liquid can also contain at its interior portion of a second phase organized with specific shapes. This is the domain of micelles, vesicles and other similar structures, all subject of intense research both of basic type and addressing practical applications.

The ordering of specific chemical systems into layers and vesicles has a counterpart in bulk liquids. A special category of liquids are called liquid crystals. They combine macroscopic properties shared by all liquids, in particular the properties of assuming the shape of the vessel in which they are put and the capability of dissolving other chemicals, with a long range order making them similar to crystals. This long range order may assume different forms, to which corresponds specific names for the liquid crystal phase, smectic, nematic, cholesteric. In general the liquid crystal are specific phases of a liquid that can also have an isotropic phase without long range order. For many substances the changes in physical parameters temperature and pressure (T, P) rule the transformation from a solid, to a liquid and then to a gaseous phase (in the order on increasing T and decreasing P), while in the substances giving origin to liquid crystals there are some additional phases between solid and isotropic liquid. The liquid crystals are a subject of intensive study, because of their quite peculiar optical and electric properties.

1 Physical approaches to the study of liquids

The distinction we introduce here between physical and chemical approaches to the study of liquid systems is mainly justified by exposition reasons. The basic principles are the same in both approaches and there is a mutual interchange of methods and procedures: both approaches grow up in harmony. Actually there is a basic difference in the motivations of the two approaches: chemistry is directly interested to the details due to differences in the chemical composition of the fluid (the study of pure liquids could be viewed in this context as an extreme case of solutions), physics arrives to consider effects due to chemical composition (and of specific solutions in particular) at the end of a longer route. The physical approaches we shall consider here do not play much attention to the chemical composition of the system (but some aspects have still to be considered).

1.1 Macroscopic approaches

Historically, and also logically, the first contribution of the physical understanding of liquids was obtained with macroscopic approaches. We shall not consider here the mechanical studies that come first in the physical enquiry about fluids, but we directly shift to the thermodynamic approach elaborated in the nineteenth century. The science of heat, thermodynamics, regards all the matter in general, and considers liquids as a specific state of aggregation of the matter, to be treated on the same footing as the others. This emphasis on the uniformity of thermodynamic laws actually is at the basis of our understanding of the phenomena of phase transformation we have quoted above, and of the transfer of models (that we shall quote in the following) elaborated for the gas phase, being simpler to study, to liquids.

Thermodynamics is a rigorous discipline, and the definition of thermodynamic functions (at the equilibrium and out of the equilibrium) must be always reminded in performing studies on liquids with theoretical tools. Even when the model is reduced to a molecular model with attention paid to the details of quantum mechanical (QM) calculations on a reduced portion of the liquid, it is wise (often necessary) to keep in mind what is the thermodynamic status of the system, and to what selection of fixed macroscopic variables (pressure, volume, temperature, energy, chemical potential) is made in assessing the models.

A second macroscopic approach of interest for us regards the electric properties of the liquid. Here again we have to go back to the nineteenth century to find the elaboration of the macroscopic theory. A large portion of liquids are poor conductors of electricity (while there are notable exceptions, the whole electrochemistry is based on the conducting properties of ionic solutions), and for this reason the dielectric behavior plays the prominent role. It is convenient to give a short summary because we shall need it later. An external electric field induces a polarization of continuum dielectric media. In the standard version of the Maxwell elaboration, the attention is focused on the dipole polarization with respect to a homogeneous electric field. A vector field, P, is defined giving the value of the dipole density; this vector is related to the two other vector fields defined by Maxwell to satisfy the two basic constitutive relations for electrostatic fields in vacuo: the electric field E and the displacement vector D (by definition E − D = 4πP.

In the simplest case (homogeneous linear dielectrics, constant electric field) the relationship P = χE with χ = (1 – ε)/4π holds. The macroscopic dielectric description of the liquids is not limited to the basic homogeneous and isotropic description we have here recalled. For anisotropic fluids (as liquid crystals) a tensorial definition of χ and of the dielectric constant ε must be introduced. There are cases in which the linear regime is not sufficient and the polarization must be described with the aid of higher order terms in the electric field. Models for specific cases call for specific modifications of the basic dielectric equation (a notable example of wide interest in chemistry is given by the Debye–Hückel model for ionic solutions). Other modifications are used for liquids spread on a metallic body.

All the examples we have given here refer to specific applications that shall be considered more in detail later, but of wider occurrence are the dynamic aspects introduced by the time dependence of the external field E(t), which gives rise to a large number of important phenomena. Relations paralleling those defined for all the static cases we have mentioned can be derived when changes in the electric field are not excessively fast (for the definition of such limits one has to make reference to a microscopic picture of matter). At a basic level, the dynamic case can be studied with the help of sinusoidally varying electric fields: E(t) = E⁰ cos ωt giving so origin to a complex dielectric response function that can be split into a real part, ε′(ω), and an imaginary part, ε″(ω). The first is a generalization of the static dielectric constant, ε′(0) = ε, and the second is called the loss factor (BÖTTCHER [1973], BÖTTCHER and BORDEWIJK [1978]). Both functions play an important role in the following of this chapter.

1.2 Microscopic approaches

The use of microscopic models based on molecules to describe liquids was introduced at the beginning of the past century. The impact has been enormous pervading all the research fields on liquids. The theory of liquid is at present a molecular theory. Two are the reasons of a so large impact: first, microscopic approaches have opened the way to the use of statistical mechanics, second, they have introduced explicit consideration of molecular interaction potentials, a subject on which chemists have a lot to say. For more than fifty years the research has been led by physicists, following the strategy to which physicists are more inclined, namely discarding the details of the molecular interaction potentials, spending more effort to establish models based on statistical mechanics and taking the model for ideal gases as starting point. Only in the late sixties of the past century aspects of more direct interest for chemistry were taken into consideration. In this chapter, we are more interested to the chemical way of describing liquids and so we shall deserve only few sentences to the physical approach, especially to introduce terms that will be used in the following chapters (among the vast literature on the argument, we quote here few titles: BALESCU [1991], REICHL [1980], MCQUARRIE [1976]).

The basic ingredient in this approach is the probability density P in the 6N-dimensional phase space of a liquid system which satisfies the invariance of a selected set of three macroscopic quantities (the number of particles N, the temperature T, the energy E, or the chemical potential μ). In the formulation of the theory of statistical thermodynamics use is made of another concept, that of ensemble of systems, an ideal collection of numerous systems, all defined by the same invariants, in close contact and exchanging among them what permitted by the invariants. These ensemble are usually indicated as microcanonical, canonical and grand canonical. We report in the table names and invariants of the systems to which are associated. In the following we shall indicate single systems with the name of the ensemble, to avoid proliferation of names. The systems that have been selected correspond, respectively, to closed and isolated, to closed but not isolated and to open systems. Several other systems may be defined, according to the experimental setting used for the study of the properties shown in Table 1.1.

Table 1.1

The time evolution of the probability density P(t) can be described in terms of a Liouville equation. The description given by the full P function is too detailed as in practice we only need the probability density to find expectation values or correlation functions for various observables. The observables generally corresponds to one- and two-body operators and so it is sufficient to use reduced densities limited at these two orders: P1 and P2, formally obtainable by integration of P on the other variables. The equations of motion for reduced densities must be however expressed in terms of a hierarchy of equations, called BBGKY, where P1 is given in terms of P2, P2 in terms of P3, and so on, until the so-called thermodynamic limit (quite large). To use the BBGKY hierarchy there is so the need of defining a closure relation, and the various physical methods differ in the choice of it. The BBGKY hierarchy was defined within the classical picture of the fluid; in quantum systems, where the Hamiltonian replaces the Liouvillian operator, use is made of the Wigner functions which lead to a similar hierarchy.

Another approach is available, based on the correlation functions. The two approaches are connected. For brevity we shall limit ourselves to stationary states of the Liouville equation, namely to systems in equilibrium. P can be expressed as a product of uncorrelated one-molecule distribution functions P1(i) modified by a correlation function gN:

We have indicated with i the set of coordinates of the ith molecule (remark that implicitly, we have here discarded the grand canonical system, that can be recovered later, and that we have not paid attention to the possible differences among the N molecules of the system). For isotropic and homogeneous fluids each P1(i) is independent on the coordinates of i, so we may express it in terms of the numeral density ρ

The correlation function is then subjected to a cluster expansion:

Let us consider the two-body correlation function g2(i, j). It can be divided into two parts, the first called direct correlation function c(i, j), describing a direct effect only depending on the interaction potential V(i, j), and the second describing many-body effects, in which i influences j through other molecules. The Ornstein–Zernicke equation, set in 1914 for spherical rigid molecules (nonadditive repulsive potential), and generalized in 1972 to rigid nonspherical molecules, is used as an approximation to connect g2(i, j) to the direct correlation function. After some manipulations (use is made of an intermediate function indicated as h(i, j)) one obtains an expansion into powers of the numeral density, in which the coefficients are given by the integration of appropriate products of direct correlation functions:

Clearly this approach presents the same problems of the equation of motion method. The same set of possible closure relations are in fact used in the two approaches. It is worth to remark that the use of these approaches are limited, with a very few exceptions regarding very simple model fluids (low pressure gases), to a calculation of g2(i, j) discarding higher terms in the cluster expansion. The same limitation is used in the computer simulations on liquids we shall present in the following of this chapter.

To close this chapter on the physical approach we consider some formal aspects of the calculation of properties. The number of properties that can be treated with physical approaches is large, but does not include many properties of chemical interest (mostly related to the behavior of specific molecules in the liquid). In parallel several properties computed with the physical approach are of little interest in chemistry. It may be worth recalling that properties object of scrutiny in the physical approach can be divided into three categories. First, the properties having a microscopic counterpart and that be defined as averages on the microscopic dynamical functions (they are often called mechanical quantities), examples are the internal energy E and the pressure. Second, the properties that have no microscopical meaning and must be defined in terms of the whole probability density (called thermal quantities); examples are the temperature T and the entropy S. Third, quantities that are fixed in value by the experiment without reference to the internal state of the system; they are called external parameters and typical examples are the number of particles N, the volume V, and the strength of external fields.

The method we are sketching here below can be applied to properties of the first and second type. We limit ourselves to consider time-independent properties in equilibrium systems. The value of a property is given as the average of the values the property has for each distribution, each multiplied by the appropriate statistical weight. The normalization factor of this sum is denoted as the partition function Z. Each kind of systems has its own partition function. There is no complete uniformity in the literature about the symbols used to indicate the corresponding partition function, here we shall use Ω, Q and Ξ for the microcanonical, canonical and grand canonical systems, respectively. The system of reference is the canonical one. The steps leading from these systems to the other will be here omitted. In classical thermodynamics the complex (and complete) set of differential relations connecting the various thermodynamical functions is expressed taking one among them as the leading term (characteristic function). In the (N, V, T) systems the characteristic function is the Helmholtz free energy A. Its expression in statistical thermodynamics is quite simple:

The most important thermodynamic functions for canonical systems are reported here below:

Among the five energy functions reported in the table the most important in chemistry is G, also called Gibbs free energy, or free energy tout court. G is the characteristic function for another type of systems: (NPT) called isobaric–isothermal, and it also plays the leading role in the definition of grand canonical systems, through a related thermodynamic function, the chemical potential (in the chemical literature the term chemical potential is often used instead of Gibbs free energy). The function A plays a limited role in chemistry, it has been introduced because the analysis in terms of the canonical ensemble is easier than using other invariants. Remark, however, that when the attention is limited to liquid systems near the standard conditions there is no difference in practice between A and G.

2 Chemical approaches to the study of liquids

The simple model of liquids used in the chemical approach considers liquids as composed by an assembly of molecules at close contact, undergoing incessant collisions. Macroscopic conditions (T and P) rule the kinetic aspects of these collisions, that are however strongly modulated by molecular interactions depending on the chemical nature of the molecules. These interactions are responsible of rare events (with respect to the number of collisions) which give origin to reactions and other phenomena of chemical interest. The same interactions are also responsible of the creation of a partial ordering of local nature, which fades away at larger distances. Chemists tend to mentally average the thermal collisions and to pay more attention to the effects of chemical interactions, both for reactions and for the partial local ordering, which affects the properties of the composing molecules, the solute in particular.

2.1 Molecular interactions

The mutual molecular interactions are of different type, some are of very short range, others exert their effects at a longer distance. In general they are nonadditive, but with some exceptions. Chemical practice has permitted to learn a lot about the effects of these forces in the various solvents and obtained some hint about their nature. A more precise classification and definition came from the study with theoretical tools of the behavior of small molecular clusters, composed by two, three or little more molecules (MARGENAU and KESTNER [1971], HOBZA and ZAHRADNIK [1988]). The classification in use for the decomposition of these interactions is reported in Table 2.1.

Table 2.1

The whole set of interactions can be described at the QM level using the usual electrostatic Hamiltonian, other terms not present in it, like spin-orbit effects or relativistic corrections play a little role in the study of liquid systems.

The complete interaction energy within a given cluster can be easily computed with standard QM tools. It is sufficient to compute the energy of the whole system, considered as a supermolecule, and then subtract the energies of the separate partners. This interaction energy can be decomposed with the help of some additional QM calculations performed with the same procedure but with deletion of some terms and/or of some steps in the solution of the resulting computational problem. The details are not reported here, suffice it to say that the analysis of such decompositions performed on small clusters at different geometries has permitted to gain a good confidence on the relative importance of these terms, and on their spatial anisotropy.

There is another approach that do not introduce the calculation of the energy of the supermolecule, it was the only approach in use before the advances in electronic computers made possible calculations on supermolecules. Use is made in this approach of the perturbation theory (ARRIGHINI [1981]). The unperturbed Hamiltonian is defined as the sum of the Hamiltonians of the isolated monomers, and the perturbation is just given by the interaction. The perturbative corrections to the energy are expressed order by order and within each order a separation in additive terms corresponding to interactions of different physical origin is performed. There are problems, however, due to the fact that the full Hamiltonian of the cluster has a higher permutational symmetry than the unperturbed one due to the exchange of electrons among partners: the wavefunction of the cluster has to be fully antisymmetrized with respect to all the electrons. This lack of antisymmetry in the exchange of electrons among partners greatly complicates the formulation of the method, and the quality of the results. Here again we do not give details about the procedures elaborated in the years to partly overcome these problems. Suffice it to say that the results obtained for lower orders of the perturbation expansion on small clusters are in fairly good agreement with those obtained with the decomposition of the supermolecule description. The basic point is that neither approaches can be directly used to get a description of the liquid around a given solute.

There are two reasons, strictly connected. The first is that experience has learned that models composed by small size clusters give a very poor description of local solvent effects (and almost nothing on the properties of the whole solution). The second is that the larger aggregates that should be used in modeling require an accurate description of the thermal motions, properly averaged to reach a consistent thermodynamic status. QM supermolecule treatments are ruled out mostly for the difficulty of getting the thermal average. QM calculations at single geometries of a large cluster are feasible with some efforts, but the spanning of the portion of interest of the conformational space, and the following determination and use of the system partition function are out of question. The device often used in the study of isolated chemical systems of looking only at minima and other topological critical points on the potential energy surface (PES) cannot be applied here. The PES presents an exceedingly number of local minima, separated by very low barriers and with flat portions of the surface, hard to treat with statistical mechanics. In practice, the concept of PES is of little utility in treating large molecular aggregates. The direct use of the explicit expressions of the various components of molecular interactions obtained either via decomposition of the supermolecule energy or perturbation theory analysis are out of question for similar reasons. The computational cost is nowadays similar to that of getting the whole supermolecule energy (in addition, the extension of such formulas to the case of many body clusters presents serious difficulties).

There is so the need of reconsidering the problem, paying attention to the physical effects produced by the various types of interactions. We shall give here a short summary of this analysis (TOMASI, MENNUCCI and CAPPELLI [2001]).

COUL. It is the interaction having the larger long range effect. It is strictly nonadditive (i.e., limited to two-body contributions) and strongly anisotropic. In most cases it is the more important contribution (for example in water solutions). It determines the most favorable orientation among partners having asymmetric charge distributions (e.g., dipoles). The formal expression for a two body interaction

suggests that simplifications of this expression can be searched with an opportune modeling of the total charge distribution ρ(r) or of the electrostatic molecular potential V (r). This last presents, as a rule for neutral molecules, positive and negative regions, and the same holds for the integrand giving origin to COUL. Acceptable simplifications must preserve the anisotropy of the property, which, as already said, is essential. One-center multipole expansions work badly if the molecule has a complex shape: better it is to pass to many center expansions. One formulation of large use in chemical computations (especially for solutions of large biological systems) consist in keeping only the first term (a point charge) of multipole expansions centered on all the atoms of the molecule. This approximation is rather grossly, and should be avoided, if possible.

IND. This term is decidedly less anisotropic than the preceding one, but it exhibits a strong nonadditivity. The integrand giving origin to IND is everywhere negative. The modeling of the interaction is generally based on multipole expansions of the electric molecular polarizability. The effect of the other molecules on the charge distribution is in general far from that of a uniform electric field assumed in the definition of the first order polarizability tensor. In spite of it, a formulation often used in computations consists in placing a single isotropic value for the polarizability placed at the center of the molecule. It is a poor approximation, partially justified for a very common solvent, water. The nonadditivity of the contribution is active also at relatively long distances: the electric field of distant molecules must be considered in the contribution at each polarizable center. In turn polarization effects, even when reduced to a single contribution, affects the total polarization. IND must be computed in an iterative way, when these models are used.

EXCH. To satisfy antisymmetry, closed electronic shells have to repel each other. This contribution is everywhere positive. It is nonadditive, but only the effect of nearby molecules play a role, being the contribution considerably short-ranged. Exchange forces can be relatively well described by a set of repulsive stiff potentials, centered on the various nuclei of the molecule, each with a spherical symmetry (typically with a R−12 decay law). The exchange forces are the interaction terms left in the oversimplified physical models, in which the molecule is replaced by a sphere, hard or with a soft repulsive potential.

DISP. This term is nonadditive, with quite moderate anisotropy. It may be formally treated as IND, making use here of dynamic multipole polarizabilities. The simple approximation used in chemical computations uses the first term in the development of the first polarizability in terms of the inverse powers of the distance (it is a R−6 term) applied again to each atom of the molecule. Dispersion forces are relatively weak and have a moderate long range effect, by far smaller that of COUL and IND. It happens however that the long range effects of the two classical electrostatic terms, coupled to their anisotropy, tend to orient molecules in the liquid in the more appropriate orientation, producing so a screening of the global effect. The screening is not active for dispersion terms, and so happens that for several of the examples of liquids we have done in the introduction, the long range interactions are almost solely ruled by dispersion forces. This effect is more evident in the presence of massive bodies.

CHTR. The effects due to the transfer of electronic charge during the molecular encounter are harder to model than the preceding ones. Of course these effects belong to the category of rare events but a large part of chemistry is based just on these rare events. In the approximated formulation of interaction effects we have outlined, electron transfer effects have not been inserted. This is a limitation in the computational procedures using this modelistic elaboration. To consider the effect of charge transfer other methods must be devised.

We have outlined a reasonable way of computing molecular interactions without passing through QM calculations. This procedure can be inserted as a basic element into other procedures performing the systematic scan of the conformational space on the desired thermodynamic ensemble, to get the properties of the liquid. This scan is given via molecular simulations, Monte Carlo (MC; according to a Gibbs picture of thermodynamic averages) or to Molecular Dynamics (MD; according to a Boltzmann picture). These topics will be considered later in this chapter. Such simulations, largely used in the last 20 years, have given a considerable wealth of information about the properties of liquids, by far more detailed and more precise than those obtainable with the physical approach. All the types of liquids we mentioned in the introduction have been at least partially examined with these approaches, and it may be said that what we know today about liquids derives from, or has been confirmed by, molecular simulations. It has been a big effort, quite rewarding. The magnitude of this effort, that has been at a good extent of methodological type, can be appreciated by looking at the final computational costs. To perform a simulation on a liquid system there is to repeat many times the computation of the geometry and of the interactions of a relatively large cluster. For this reason people is compelled to use descriptions of the molecular interactions as simple as possible. In addition, this approach does not give the essential information for which chemists undertook studies on liquids: the effect of a solvent on a chemical reaction and on the molecular properties of the solute. To have them a QM description at high level is in fact necessary.

An approach alternative to the simulation on simplified semiclassical interaction potentials exists, and it is able to reach the requested high quality in the description of the molecular units of interest (TOMASI and PERSICO [1994]). This approach is the main subject of this chapter. It is still based on the analysis of the various interaction terms, but in a different way. The basic consideration is that the thermal average performed as the final step in the above described methods can be avoided by replacing the molecular discreteness of the molecular distribution with a continuous distribution function. It is possible to formally define continuous solvent response functions corresponding to the various interaction terms, and satisfying the macroscopic