Methods for Determining Partial Charges

Question

I want to run classical molecular dynamics simulations of a periodically replicated surface (rutile $\ce{TiO2}$ with grooves). In order to do so, I first need to solve for the partial charges residing on each atom of the surface. What electronic structure methods are available for determining the partial charges on a surface?

I am vaguely familiar with Mulliken population analysis, but my understanding is that this method is not the most accurate. What other methods are commonly used, and how accurate are each of these methods?

Finally, what are some good (and free) packages for performing such calculations? I've used SIESTA in the past, is that well-suited for determining partial charges of a surface?

Answer 1

There are lots, and I mean lots of methods and software programs to produce partial charges. See Wikipedia for a small (incomplete list)

Let's start with the basics. The idea of a partial atomic charge, while useful for concept, can not be defined uniquely. Quantum chemical methods (whether wavefunction or DFT) produce some sort of electron density. Different schemes divide up that electron density in different ways.

I like Cramer and Truhlar's categories:

• Class I: Some sort of intuitive or empirical approach (i.e., non-QM). These include methods based on experimental dipole moments and atomic electronegativities, including Gasteiger-Marsili partial charges and other "electronegativity equalization methods" (EEM).
• Class II: Partitioning using wave functions / orbital schemes (including Mulliken charges while very efficient to compute, are basis set dependent.
• Class III: Partitioning the electron density or electrostatic potential (e.g., Hirshfeld, density-fitting schemes, etc.
• Class IV: Semiempirical mapping from Class II or Class III (preferred) to match experimental data like dipole moments.

Almost any modern quantum package will include several schemes to fit partial charges. My personal preference is for something like Hirshfeld charges (i.e., from the electron density, so not basis-set dependent), or electrostatic potential fitting schemes like Merz-Kollman or CHelpG.

I think the latter (electrostatic potential fitting) will be better generally for molecular dynamics, since you're attempting to produce point charges that will best represent the quantum mechanical electrostatic potential.

Answer 2

Kinda obvious, but still.

First, you have to decide, what face you are working with. $\ce{TiO2}$ may exhibit several different faces, you'll have to consult literature to find out which are the most frequent.

Second, you have to consult literature to find which surface species exist on the surface in your conditions. Generally, oxides tend to adsorb water forming hydroxide groups and some other small molecules.

After that you have to build a model slab of the surface. Generally, two-three layers of atoms close to surface are 'relaxed' and two-three 'deeper' ones are fixed in positions, modelling unperturbed crystal structure. However, ideally you should perform several computations, gradually increasing amount of layers in slab until the characteristic you considering stabilize.

It is generally recomended to employ periodic approach using plane-wave sets, possibly in DFT method, as it greately reduces memory requirements and amount of atoms (and electrons) considered. Given you care to use the charges to model electrostatic potential, Mulliken and Natural charges go to garbage. Ideally you should use full-scale electrostatic potential. However, if it is not an option, analog to Gaussian Pop=CHelpG should be used, that fit atomic charges to produce electrostatic potential closest to one produced for given electronic distribution. Not aware of other programs implementing similar schemes, so you have to find one on your own.