Transferability of Atomic Multipoles in Amino Acids and Peptides for

ARTICLE pubs.acs.org/JPCA

Transferability of Atomic Multipoles in Amino Acids and Peptides for Various Density Partitions Magdalena Woi
nska and Paulina M. Dominiak* Department of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warsaw, Poland

bS Supporting Information ABSTRACT: Multipole expansion of electron density distribution is an efficient tool for evaluating the energy of interactions in molecules. In as much as atoms in macromolecules such as peptides are modeled with certain types of atoms derived from small organic molecules, investigating transferability of atomic multipoles for various partitions of molecular electron density is an important issue. In this study, multipole moments up to hexadecapoles for types of atoms present in selected amino acids, as well as di- and tripeptides composed of these amino acids, are computed using three density partitions: HansenCoppens aspherical pseudoatoms formalism, Hirshfeld’s stockholder partition, and Bader’s atoms in molecules theory. Electron density of relevant molecules is derived in a procedure including molecular wave function ab initio calculation for isolated molecules in geometry from X-ray measurements, calculation of theoretical structure factors for molecules put in a pseudocubic cell, and multipole refinement as in crystallographic data processing and computation of multipoles. The results were compared to calculations of multipole moments in AIM and in stockholder density partitions obtained directly from molecular wave functions. The presented comparison does not point unambiguously to any particular influence of multipole refinement on moments obtained from these two partitions. The advantage of stockholder partitioning in terms of transferability of atomic multipoles is affirmed. AIM and pseudoatoms provide slightly less transferable multipoles of lower order. Higher rank of multipole expansion reveals a transferability improvement in the case of AIM and meaningful deterioration for pseudoatoms.

’ INTRODUCTION A method commonly applied in molecular mechanics for approximating the potential energy of large molecules is force fields. Many of them restrict the electrostatic component of the energy only to the interactions of atomic charges centered on the nuclei. Often such an approach is too simplified; therefore, in new generation force fields,1
3 higher multipole moments are also included. Creating chemical databases of atomic charge densities based on highly transferable multipole moments4
7 allows to reach a good model of electrostatic properties of complex molecules at low computational cost. This method can also be employed in crystal structure refinement.8
10 Electron Density Partitions. The possibility of assigning to atoms in a molecule additive properties, such as multipole moments, which summed up will reproduce corresponding properties of the whole molecule, evokes the conception of transferability. Atoms in a large molecule may be divided into some types with respect to certain parameters characterizing them. To simplify the calculations, properties of all atoms of the same type may be considered equal and obtained from a suitable database. So that a complex molecule could be modeled properly, atom types must be defined in the way providing the highest similarity of representative atoms, which means that chemical qualities must be highly transferable among atoms belonging to the same type. r 2011 American Chemical Society

Because properties of atoms may be introduced as functionals11 of atomic electron density, their transferability is dependent on how molecular density is shared by individual atoms, that is, on the method of density partition. Herein, three commonly used density partitions are considered. Two of them are fuzzy, which means that densities of atoms overlap in space, unlike the third one which is discrete. A brief description of each one is given in this section. The aspherical pseudoatoms formalism is a fuzzy partition based on the multipole model of electron density implemented by Hansen and Coppens.12 In this method, molecular electron density is a sum of atomic contributions called pseudoatoms. The density of the i-th pseudoatom is given by the expression r Þ ¼ Fi, c ð B r Þ þ Pi, v Fi, v ðki B r Þ þ Fi, d ðk0i B rÞ Fi ð B

ð1Þ

where Fi,c and Fi,v are spherically averaged Hartree
Fock core and valence electron densities13 and Fi,d is an aspherical deformation of valence density. Parameter Pi,v denotes population of the spherical valence shell, and ki and k0 i are screening constants. The first two terms of eq 1 stand for the Received: April 29, 2011 Revised: August 26, 2011 Published: September 26, 2011 1535

dx.doi.org/10.1021/jp204010v | J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

spherical part of atomic electron density, while the last term, r ), represents aspherical deformation of valence denFi,d(k0 i B sity and is equal to
 r r Þ ¼ Σl Rl ðk0i rÞΣml¼
m Plm, i ylm B Fi, d ðk0i B ð2Þ r Rl are Slater-type radial functions, Plm,i are the multipole populations and ylm denote real spherical harmonics further defined by eq 9. This model is frequently used for multipole refinement of high resolution X-ray crystallographic data. Core electron density Fi,c(rB) is always frozen, whereas the population parameters Pi,v and Plm,i, as well as the screening constants ki and k0 i are adjusted in the refinement procedure. The stockholder partition introduced by Hirshfeld14 defines atomic electron density utilizing the notion of a reference state. The electron density of a reference state is the spherically averaged ground-state electron density of an isolated atom. The sum of the densities of atomic reference states centered on the respective nuclei forms the density of a promolecule. In the stockholder partition the participation of the i-th atom Fi(rB) in molecular electron density F(rB) in the point of space denoted by vector B r is proportional to the corresponding reference state’s Fi,ref(rB) contribution to the promolecule’s density Fprom(rB) and can be expressed as rÞ ¼ Fi ð B

Fi, ref ð B rÞ Fð r Þ Fprom ð B rÞ B

Fprom ð B r Þ ¼ Σi Fi, ref ð B rÞ

ð3Þ

ð4Þ

Hirshfeld’s model can also be applied for crystallographic analysis. The third method of density distribution results from Bader’s theory of atoms in molecules (AIM).15
17 This method is a discrete partitioning in which the electron density of each atom is connected with a separate region of space, called the atomic basin (denoted Ωi). This area is defined as the volume limited by the so-called zero-flux surface: r ∈ R3 : ∇Fð B r Þ3 B nð B r Þ ¼ 0g ∂Ωi :¼ f B

into a series called a multipole expansion:18,19 1 Vð B r Þ ¼ Σ∞ Σl ε l¼0 m¼

where ε is the dielectric constant and B r is a vector determining the position of a point at which the potential is calculated. If the point designated by vector B r is localized far from the area of nonzero electric charge (r . r 0 ), the potential can be developed

qlm Ylm ðθ, jÞ 2l þ 1 r l þ 1

ð7Þ

where θ and j are angular coordinates of vector B r and Ylm(θ,j) are spherical harmonics. Spherical multipole moments qlm are given by the formula Z  dB r 0 Fel ð B r 0 Þr 0l Ylm ðθ0 , j0 Þ ð8Þ qlm ¼ The complex conjugate can be skipped if the basis of real spherical harmonics is used. Real spherical harmonics ylm are created as linear combinations of complex spherical harmonics Ylm: 9 8 > > Y if m ¼ 0 lm > > > > > = < p1ffiffiffi ðY þ ð
1Þm Y Þ if m > 0 > lm l;
m ð9Þ ylm ¼ 2 > > > > 1 m > > > ; : pffiffiffi ðYl;
m
ð
1Þ Ylm Þ if m < 0 > i 2 Real spherical harmonics with m > 0 are harmonics of the cosine type (denoted also yl|m|c), and those with m < 0 are harmonics of the sine type (denoted yl|m|s). This notation is used for the components of multipoles presented in this work. Multipole expansion allows to obtain fine results for calculations of interaction energy if the series is truncated at an appropriate moment. Atomic multipole moments can be used as properties characterizing atom types present in a certain class of chemical compounds and their values can be collected in a databank of atoms.20
22 Structure Factor. X-ray diffraction is a very useful tool for experimental determination of molecular electron density. The intensity of X-radiation scattered by electrons in a crystal is proportional to a squared magnitude of a static structure factor (or structure factor for short). Structure factor F(h B) is a Fourier transform of electron density, therefore the electron density can be calculated on the basis of a structure factor as a Fourier summation: Fð B rÞ ¼

ð5Þ

where the surface bounding the atomic basin of the i-th atom (Ωi) is symbolized by ∂Ωi. rF(rB) is the gradient of the electron density at the point B r and B n (rB) (called the normal vector) is a r . The vector perpendicular to the surface ∂Ωi at the point B electron density gradient at each point of the zero-flux surface is perpendicular to the normal vector at this point (eq 5). It means r and, therefore, the that rF(rB) is tangent to ∂Ωi at each point B flux of rF(rB) through the whole zero-flux surface vanishes. Multipole Moments. Electrostatic potential is an important property in the description of interactions in a crystal. The component of potential resulting from the electronic charge density Fel(rB) is given by Coulomb’s law:18 Z Fel ð B r 0Þ 1 Vð B rÞ ¼ ð6Þ dr0 4πε j B r
B r 0j B


l

1 Σ Fð hBÞe
2πi Bh Br V Bh

ð10Þ

where B h is a vector from the reciprocal space, describing an X-ray reflection obtained as a result of scattering on a crystal lattice; B r is a vector from the real space describing position in a crystal and V is the volume of a unit cell. The static structure factor, as opposed to the dynamic structure factor, does not describe thermal fluctuations of electron density.

’ COMPUTATIONAL DETAILS Selected Compounds. A total of 10 types of amino acids were chosen for the study (alanine, serine, valine, threonine, isoleucine, leucine, aspartic acid, asparagine, phenylalanine, and tyrosine) together with some di- and tripeptides composed of these amino acids. Crystallographic structures of all the compounds were retrieved from the Cambridge Structural Database.23 Only data not containing disorder or metallic ions were selected, structures with an R-factor higher than 5% were also avoided. Only molecules in zwitterionic form are considered. Due to the demand for a statistically meaningful number of atoms representing each type, additional data from crystallographic structures of small organic molecules was used. 1536

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Table 1. Distinguished Carbon Atom Types (in Red) with Their Local Reference Framesa

a

Substitution of the nearest neighbors is given in brackets. Number of atoms representing each type is given for calculations with crystallographic refinement (density) and for AIM calculations from wave function (AIM, wfn). A line connecting two atoms may denote either single or double bond.

Theoretical Calculations. Molecular geometries resulting from an X-ray experiment were used. X
H distances were extended to their standard values from neutron measurements,24 yielding more accurate positions of hydrogen atoms. Ab initio calculations for isolated molecules were performed with program Gaussian 0325 using density functional theory (DFT) with application of Becke’s three-parameter hybrid method26 with the nonlocal correlation functional of Lee, Yang, and Parr (B3LYP).27 The standard split-valence double-exponential 6-31G** basis set including polarization functions28 was used. Complex static valence-only (calculated on the basis of valence electron density) structure factors in the range of 0 < sin θ/λ < 1.1 Å
1 were obtained by the analytic Fourier transform of the molecular charge densities for reciprocal lattice points corresponding to a pseudocubic cell with 30 Å long edges.29 Aspherical Pseudoatoms Refinement. The data were refined with the XD2006 program suite.30 Model structure factors calculated from aspherical pseudoatoms in the Hansen-Coppens formalism12 were fitted to the theoretically computed structure factors. Phases were fixed at the theoretical values. Radial screening factors k and k0 were refined independently for every atom, except chemically equivalent hydrogen atoms, for which both of these coefficients were equal. The value of the total charge of a

molecule was fixed and no local symmetry constraints were imposed. The refinement procedure included multipoles up to hexadecapole level (l = 4) for all of the atoms apart from hydrogen. In the latter case, the expansion series was truncated at quadrupoles (l = 2) and only multipoles with m = 0 (bonddirected) were used. Radial functions in the pseudoatom formalism were taken as simple Slater functions with the coefficient in the exponent deduced from single-ζ wave functions.31,13,32 Subsequently, multipole moments for three described methods of electron density partition were calculated. Computations were performed with module XDPROP of the XD program. Atomic Cartesian unabridged multipole moments33,34 in global reference frame were calculated for atoms derived from aspherical pseudoatom partition,12 stockholder partition14 and via integration over atomic basins15 (Cartesian unabridged multipole moments are defined by the following equation: Z 1 d3 rFð B r Þrα1 rα2 :::rαl μa1 α2 :::α1 ¼ l! in which rαi stands for one of the Cartesian coordinates x, y, or z (i = 1, 2, ..., l), the coefficients α1α2 ... αl denote the component of the multipole, and the integer l describes the order of the 1537

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Table 2. Distinguished Hydrogen, Nitrogen, and Oxygen Atom Types (in Red) with Their Local Reference Framesa

a

Substitution of the nearest neighbors is given in brackets. Number of atoms representing each type is given for calculations with crystallographic refinement (density) and for AIM calculations from wave function (AIM, wfn). A line connecting two atoms may denote either single or double bond.

multipole). Then, the multipoles were transformed to the local spherical coordinate systems of relevant atoms. In the case of multipole moments computed from aspherical pseudoatoms, transformation of the coordinate system was done by a locally modified version of XDPROP. Such a procedure was performed to follow the scheme of crystallographic data processing commonly applied in experimental charge density analysis. To test the influence of pseudoatom refinement on the multipoles, calculation of AIM moments directly from the Gaussian wave function was performed with program MORPHY98.35 Additionally, stockholder moments up to quadrupoles were calculated straight from the wave function using program TONTO.36 In each case of the computational method, components of multipoles were averaged over atoms representing a given type. Standard deviation resulting from the averaging was accepted as a measure of transferability level. Atom Type Definition. In the amino acids and peptides under research, 27 types of atoms were distinguished. Criteria leading to such division are similar to those used in databanks of transferable pseudoatoms,5,37,9 that is, (1) chemical element type, (2) number and type of the nearest neighboring atoms, and (3) substitution of the nearest neighbors. Aromaticity is also included in atom type definition. Chirality is not taken into account, as it only influences the choice of local reference frame of an atom

(right-handed for configuration S and left-handed for R). All the types of atoms, along with orientations of two axes of the reference frames in which are expressed atomic multipoles, are given in Tables 1 and 2. The third axis is perpendicular to the other two and is oriented to create a right-handed frame (with the already mentioned exception of chirality). For hydrogen atom types, the axis y is aligned in the plane defined by the hydrogen atom, its nearest neighbor, and the next nearest neighbor with the highest atomic number and the smallest substitution. Local symmetry groups and numbers of atoms belonging to each type are juxtaposed in Table S1 in the Supporting Information.

’ RESULTS AND DISCUSSION In this section, analysis of the calculated atomic multipoles is undertaken to investigate which of the three considered density partitions provides the highest level of transferability. Evaluation of Reliability of the Results. Utilized criteria of credibility of the refinement results are based on R-factor of structures after refinement, statistical parameter Mean(Shift/su)30 (an average over all refined parameters for value of least-squares shift, divided by standard uncertainty), and values of residual density (calculated from the difference between structure factor obtained from the model and theoretical structure factor). 1538

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Values of monopole moment M with standard deviations for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement). The values corresponding to individual atom types are spread apart for better transparency.

R-Factor for various structures oscillates between 2.5 and 3.9%, which is a fine result. It is accepted that a satisfactory value of the parameter Mean(Shift/su) should be less than 0.001, which is the case for all the crystal structures. Residual density values are limited to the range of
0.15
0.1 e/Å3. Positive values are located in the positions of free electron pairs of oxygen and nitrogen atoms and some bonds. Negative values of residual density appear in the positions of nuclei (which is characteristic of Hansen-Coppens multipole model), and sometimes negative values up to
0.1 e/Å3 appear in the positions of bonds linking carbon atoms with heteroatoms. For atomic basins partition, multipoles are obtained through integration over atomic basins. Correctness of integration is evaluated on the basis of summed electric charge and volume of atomic basins of all the atoms in a molecule and the value of integrated Lagrangian. Charges of atoms sum to the total molecular charge for each of the investigated molecules. Volumes of atomic basins in a molecule should sum to the volume of a unit cell. Unfortunately, this condition is not satisfied, as the unit cell applied in the calculations is very huge (with 30 Å long edge). Module of integrated Lagrangian should not be higher than 1.0  10
4 a.u. for nonhydrogen atoms and 1.0  10
3 a.u. for hydrogen atoms. This requirement is satisfied for hydrogen atoms, but in case of nonhydrogen atoms, the permissible value is often exceeded (it sometimes amounts up to 0.01 a.u.). It may lead to some inaccuracy of the results. In the case of AIM moments calculated from the wave function, the quality of integration is even worse. There are quite a few hydrogen and oxygen atoms for which the integrated Lagrangian reaches 0.01 a.u. For carbon and nitrogen atoms, the upper limit of integrated Lagrangian is usually exceeded and it may attain 1 a.u. for carbon and 0.1 a.u. for nitrogen.

The summed electric charge of atoms and the summed volume of atomic basins cannot be checked, as in the case of more complicated interatomic surfaces integration could not proceed. Reliability of stockholder moments calculated from wave function can be estimated on the basis of summed electronic charge of atoms, which in this case sums very accurately to the total number of electrons in a molecule. Analysis of Multipole Moments. Multipole moments obtained with various methods of electron density partitioning and expressed in suitable local reference frames are compared in terms of their values (obtained via averaging over the results for atoms representing the same type), standard deviations of values (which are the measure of transferability), and compliance with constraints imposed by local symmetry of individual atoms. The last criterion allows to determine whether a particular density partition reveals lower symmetry than originally assumed and whether there are any differences between partition methods in this regard. AIM and stockholder moments computed with a passage through the multipole refinement are also compared to the corresponding moments calculated straight from the wave function in a suitable density partition. Although the numbers of atoms representing particular types are substantially different from one another, in the course of analysis, this fact turns out not to have any particular influence on the standard deviations of the multipoles. Monopole Moment. The most outstanding result is that stockholder partitioning in case of the monopoles is a method guaranteeing the highest level of transferability. Aspherical pseudoatoms partition yields slightly worse transferability than the atomic basins method, which in turn, as the only discrete partition, produces values of monopoles that are often meaningfully different from the values obtained with the two other partitions. The 1539

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Figure 2. Values of standard deviations of monopole moment M for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement).

Figure 3. Values of dipole component D10 with standard deviations for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement). The values corresponding to individual atom types are spread apart for better transparency.

group of atom types for which the AIM monopoles are the least transferable includes mostly carbon atoms connected with two electronegative atoms (oxygen, nitrogen) or being members of an aromatic ring. The values of monopoles for individual atom types and their standard deviations are given in Figures 1 and 2.

Dipole Moment. In this case, stockholder partitioning still provides the best transferability. The two other methods give alternating results, being a slight disadvantage to the AIM moments. An exceptional group is the hydrogen atom types, for which the dipole component D10 (the only one used in crystal 1540

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Figure 4. Values of standard deviations of dipole component D10 for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement).

refinement) has the highest values of standard deviation in stockholder partition, whereas the two other methods return comparable results. Components D11c and D11s for hydrogen atom types are the least transferable in stockholder partition and perfectly transferable if calculated from pseudoatoms, as this partition is favored by the procedure of multipole refinement (setting D11c and D11s to zero). Among other groups of atoms, some additional tendencies are visible. In the case of carbon atom types, the three density distributions may be arranged in the following order (from the best to the worst transferability): stockholder, pseudoatoms, atomic basins. For oxygen atom types, a suitable series would be the following: pseudoatoms, stockholder, atomic basins. The results are the most variable in the case of nitrogen atom types, for which stockholder partition mostly produces the lowest standard deviations, and AIM in the majority of examples yields the worst transferability. As for the values of components of the dipole moment, the atomic basins partition constantly manifests the most significant differences in comparison with other methods. For nonzero moments, it often yields values of opposite sign than the two other methods. Values and standard deviations for component D10 are presented in Figures 3 and 4. Other components can be found in the Supporting Information. Quadrupole Moment. For the quadrupole moments there seems to be only a small change in the proportions of standard deviations for the pseudoatoms and atomic basins partitions. Transferability in the case of aspherical pseudoatoms deteriorates for some carbon atom types. In an overall view, the stockholder method continuously results in the best transferability, although for various groups of atoms, proportions of standard deviations prove the advantage of different methods. For hydrogen atom types, the symmetry-forbidden quadrupole components, as for dipoles, are characterized by perfect transferability if obtained from pseudoatoms and are the least transferable when stockholder partition is employed. Similarly, Q20, the only component used in multipole refinement, is characterized by an increase in

standard deviations in the same sequence as it is for the D10 component of the dipole moment (atomic basins < stockholder < aspherical pseudoatoms). For oxygen atom types, this order is reversed for all the quadrupole moment’s components (as it is for the dipoles). For carbon atom types, the results speak in favor of stockholder method and are alternating for the two other partitions. In the case of nitrogen atom types, the worst results are returned by atomic basins and it is difficult to conclude which of the remaining methods is the better one. For none of the partitions are the mean values of the quadrupoles distinguished in a particular way. Values and standard deviations for the component Q20 are presented in Figures 5 and 6. Other components can be found in the Supporting Information. Octupole Moment. In this rank of multipole expansion, a new trend can be observed; the level of transferability in the pseudoatoms density distribution noticeably worsens for a meaningful number of atoms. This effect is the most outstanding in the group of carbon atom types, in which standard deviations for the remaining two methods are much lower (with a small advantage of stockholder partition). In the case of nitrogen atom types, pseudoatoms compared to AIM is generally a less favorable partition (particularly for atom type 23, four-valent nitrogen) and stockholder constantly is the most advantageous method. On the other hand, the specified changes do not apply to oxygen atom types, for which the sequence pseudoatoms < stockholder < atomic basins still describes a decreasing level of transferability. Also, in the case of hydrogen atom types, the superiority of pseudoatoms imposed by the multipole model of density persists, but the advantage of atomic basins over stockholder partition is less pronounced. For nonzero components of octupoles, those resulting from the AIM method have opposite signs than those obtained from the other partitions. At the same time, nonzero components for pseudoatoms have the most distinctive values and the highest standard deviations. However, transferability of pseudoatoms for zeroing octupoles is as good as for the 1541

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Figure 5. Values of quadrupole component Q20 with standard deviations for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement). The values corresponding to individual atom types are spread apart for better transparency.

Figure 6. Values of standard deviations of quadrupole component Q20 for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement).

remaining methods, the influence of the chemical environment on nonzero octupoles is strongly underlined for pseudoatoms. Values and standard deviations for the component O30 are given in Figures 7 and 8. Other components can be found in the Supporting Information.

Hexadecapole Moment. The trends present in transferability of the octupoles persist. As a general rule, pseudoatoms manifest very high standard deviations in comparison with the other partitions. Variations in transferability between the stockholder and atomic basins methods become less pronounced. For nonaromatic 1542

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Figure 7. Values of octupole component O30 with standard deviations for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement). The values corresponding to individual atom types are spread apart for better transparency.

Figure 8. Values of standard deviations of octupole component O30 for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement).

carbon atom types, both methods give similar results. For aromatic carbon atom types, stockholder partitioning yields better transferability for the majority of components. Hexadecapoles in pseudoatoms partition have very poor transferability for carbon atom types. The situation is not so dramatic for nitrogen atom types, where the pseudoatoms partition is not much more unfavorable than AIM (again with the exception of four-valent nitrogen nr 23). Finally, for oxygen atom types transferability still worsens in the sequence pseudoatoms < stockholder < atomic basins. Continuously, the respective sequence for hydrogen atom types is pseudoatoms < atomic basins < stockholder. For the last two groups of atomic types the differences in standard deviations between individual partitions are very small.

As for the actual values of the hexadecapolar components, for pseudoatoms many of them are extraordinarily high and substantially different from the results obtained with other methods. Still it can be observed that for some nonzero components atomic basins yield opposite values than the remaining partitions. Values and standard deviations for the component H40 are given in Figures 9 and 10. Other components can be found in the Supporting Information. Analysis of Local Symmetry. Among the distinguished types of atoms, representatives of five symmetry point groups have been identified. These are 1, m, mm2, 3m, and ∞ (cylindrical symmetry, used for hydrogen atoms). Symmetry elements belonging to each of these groups together with some examples can 1543

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

Figure 9. Values of hexadecapole component H40 with standard deviations for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement). The values corresponding to individual atom types are spread apart for better transparency.

Figure 10. Values of standard deviations of hexadecapole component H40 for distinguished types of atoms in the analyzed density partitions (wfn, calculations straight from wave function; density, calculations after multipole refinement).

be found in Table S2 in the Supporting Information. The highest possible local symmetry was preliminarily assumed for each atom. If an atom is characterized by a certain point group, the electron density of this atom has symmetry corresponding with this group. Therefore, atomic multipoles calculated on the basis of this atom’s electron density also have to comply with the restrictions resulting from the same symmetry group. These restrictions impose the requirement that those components of multipoles (further called “symmetry-forbidden”), which would violate local symmetry, must vanish. This can be the observed only if the local reference frame of an atom is properly oriented (one axis is aligned along the symmetry axis of an atom and, as far as possible, all the axes are either perpendicular or parallel to the

symmetry planes of an atom), as values of multipoles are dependent on the choice of the reference frame. For the orientation of the axes, as presented in Tables 1 and 2, the following components of multipole moments must equal zero for the specified symmetry groups: 1 — none; m - D10, Q21c, Q21s, O30, O32c, O32s, H41c, H41s, H43c, H43s; mm2 - D11c, D11s, Q21c, Q21s, Q22s, O31c, O31s, O32s, O33c, O33s, H41c, H41s, H42s, H43c, H43s, H44s; 3m - D11c, D11s, Q21c, Q21s, Q22c, Q22s, O31c, O31s, O32c, O32s, O33s, H41c, H41s, H42c, H42s, H43s, H44c, H44s; ∞ - all the components except M, D10, and Q20. It can be concluded that for nonhydrogen atom types, all three considered methods of density partition comply with the presented symmetry constraints quite precisely. Most of the 1544

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A specified symmetry-forbidden components in each partition amount to zero within one standard deviation of the mean; for some of them, this range must be extended to two standard deviations of the mean. Such results suggest that the initial choice of local reference systems was proper, as well as that all the investigated partitions manifest similar ability to reproduce the symmetry of atomic electron density. As for hydrogen atom types, pseudoatoms are a privileged partition, because none of the symmetry-forbidden multipole components is included in multipole refinement. In fact, those components in the pseudoatoms partition calculated on the basis of the refined form factors are precisely equal to zero with achievable numerical accuracy. The two other partitions usually comply with the restrictions, too; nevertheless, in the case of certain atom types, there appears some problematic components that are equal to zero only within the range of four or five standard deviations of the mean. More precisely, there are four such cases for the atomic basins distribution and six for the stockholder, which still is not particularly alarming. It was proven5 that the electron density of the carbonyl oxygen (type 24) in amino acids is deformed due to the interaction of its free electron pairs with other oxygen’s or nitrogen’s electrons, which results in symmetry reduction. For this atom type, in the case of stockholder partition, the symmetry-forbidden component D11c was equal to zero only in the extended range of three standard deviations of the mean. However, this individual result does not allow to deduce that in this study an effect of symmetry perturbation for carbonyl oxygen was observed. AIM Moments Directly from Wave Function. To verify the accuracy of computing multipole moments after the procedure of multipole refinement, multipoles (up to quadrupoles) were also calculated directly from the wave function, without employing the formalism utilized in crystallography. The AIM density partition was used in calculations. For some atoms present in the analyzed molecules, it was impossible to conduct successful integration over the atomic basin for the reason of numerical difficulties. Therefore, numbers of atoms representing certain types decreased (the numbers are given in the Supporting Information). Fortunately, this decline concerned mainly atom types with numerous representatives. Most components of multipoles calculated on the basis of the refined density do not differ significantly from the values obtained straight from the wave function. Symmetry constraints are satisfied with equal precision for both methods and the frequency of symmetry violations by multipoles of hydrogen atom types is also comparable. However, standard deviations of multipole components from the latter computational method are much higher in the case of carbon atom types and nitrogen nr 21. This may suggest that passing through the multipole model reduces the variability of atomic properties resulting from chemical environment of individual atoms. As previously stated, in the maps of residual density some minima and maxima appear in positions of the nuclei, certain chemical bonds, and free electron pairs. It may happen that even apparently small values of residual density cause considerable differences between multipoles calculated from refined and real electron density. The second factor causing high standard deviations may be problems with integration of properties over the atomic basin. These problems are visible in crash of integration in the case of certain atoms and huge values of Lagrangian integrated over the atomic basin (much higher than for calculations based on electron density). Therefore, it is impossible to conclude how the multipole model changes calculated properties without suitable investigations for stockholder moments.

ARTICLE

Stockholder Moments Directly from the Wave Function. Additional calculations of stockholder moments (up to quadrupoles) from the wave function were performed to verify the results from crystallographic refinement. In this partition, standard deviations of individual multipole components computed with and without a passage through multipole refinement are much the same in the majority of cases. A similar effect is observed for the mean values of multipole components. Moreover, constraints imposed on the multipoles by local symmetry of atom types are fulfilled very precisely. This suggests that loss of information about electron density due to applying the multipole model is not as dramatic as would result from the analysis of AIM moments and that numerical inaccuracy of integration cannot be disregarded. Summary of the Results. For each of three considered density partitions, it is possible to distinguish certain groups of atom types, for which particular partition produces the most or the least transferable multipoles. These groups are - Nonaromatic carbon atom types: multipoles with the highest level of transferability are yielded by the stockholder partition in the case of the lower-order moments, but for octupoles, AIM becomes a competitive partition, and finally, for hexadecapoles, this method becomes slightly more preferable, but still the differences are very small. The worst transferability is obtained for pseudoatoms (except dipoles and some quadrupoles), especially for higher-order multipoles, standard deviations are extremely huge. - Aromatic carbon atom types: pseudoatoms are generally the most unfavorable distribution, except for monopoles and dipoles, for which atomic basins are the inferior method. The stockholder partition produces the lowest standard deviations for the majority of atomic types belonging to this group. - Hydrogen atom types: for the only components taken into account in multipole refinement (M, D10, and Q20) utilizing pseudoatoms consequently causes the poorest transferablility. As far as the two other methods are concerned, the atomic basins partition is the best for hydrogen atoms linked with heteroatoms and stockholder is the most reliable for other hydrogen atoms. As for the remaining multipole components, pseudoatoms unsurprisingly provide excellent transferability; however, in contrast, stockholder distribution performs slightly worse than atomic basins. - Nitrogen atom types: the stockholder distribution is undoubtedly the most advantageous method. As for the two remaining ones, the results are alternating. On the one hand, pseudoatoms result in the highest standard deviations for monopoles; on the other hand, atomic basins cause transferability deterioration of dipoles and quadrupoles. The results for octupoles and hexadecapoles reveal a certain advantage of atomic basins, though standard deviations for pseudoatoms are not as large as in the case of carbon atom types. - Oxygen atom types: this is the most distinctive group of atomic types and the only one in which apsherical pseudoatoms invariably provide the best transferability. Atomic basins produce in turn the least transferable multipoles. Only standard deviations of monopoles grow in the order stockholder < atomic basins < pseudoatoms.

’ CONCLUSIONS In this work, transferability of atomic multipole moments derived from three methods of density partition (pseudoatoms, stockholder, and AIM) was compared. Multipole moments were 1545

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A calculated on the basis of theoretical structure factors subjected to multipole refinement for molecules in experimental geometries. Such a computational procedure was designed to follow the approach commonly employed in analysis of diffraction data. It is also compatible with the idea of creating a databank of atomic charge densities, a useful tool in charge density crystallography. Atoms present in certain amino acids and small peptides were chosen for the investigations in view of application for reconstructing properties of biological molecules. In particular, atomic multipole moments of high transferability are essential in the generation of new force fields, which include higher terms of multipole expansion in calculations of electrostatic interaction energy. In the selected molecules, particular types of atoms were distinguished as types, including atoms, of presumably similar properties. The accepted measure of transferability of the multipoles of individual atom types was standard deviation obtained from averaging individual components within all the representatives of a given type. The main conclusion is that the method guaranteeing the highest transferability of multipoles for the majority of atoms is the stockholder partition. This result is not surprising if one takes into account the fact that the two remaining partitions have features that may influence negatively their ability to describe atomic electron density. AIM is a discrete partition that limits atomic density to some fragment of space and that is vulnerable to numerical problems with integration over the atomic basin. Although usually the level of transferability achievable in nonoverlapping partitions is worse than in the fuzzy ones, the AIM method yields multipoles that are only slightly less transferable than the stockholder moments. Pseudoatoms, in turn, have problems with localization38,39 and, consequently, too much of electron density far from the nucleus may be attached to an atom or a meaningful part of it close to the nucleus may be missing. This effect, as already stated, is visible as nonzero values of residual density, which may significantly change multipoles. Lower order moments computed from pseudoatoms have standard deviations not varying substantially from those obtained with the other methods. Considerable decrease of transferability achieved with pseudoatoms is observed for octupoles and hexadecapoles. One can suspect that computing moments on the basis of electron density obtained from structure factors subjected to multipole refinement causes certain loss of information and can make the pseudoatoms partition privileged in comparison with the others. It is true for those components of multipoles which in the case of hydrogen atom types are set to zero in the refinement procedure. However, no similar conclusion can be drawn for other components of multipole moments. Furthermore, comparing calculations that run directly from the wave function with the results obtained after multipole refinement does not point to any significant differences in the symmetry of atomic multipoles computed on the basis of these two methods. The most important result of the research is that the stockholder partition provides the highest transferability level of atomic multipole moments with reference to the investigated amino acids and peptides. AIM returns a little worse results, but is computationally more demanding. Finally, the aspherical pseudoatoms partition, which is in fact a fit to electron density, is less reliable for higher terms of multipole expansion.

’ ASSOCIATED CONTENT

bS

Supporting Information. Plots of the remaining components of multipole moments and their standard deviations;

ARTICLE

a table containing numbers of atoms representing each type in various computational procedures; a table with a description of point groups mentioned in the paper. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT Support of this work by the Polish Ministry of Science and Higher Education through Grant Iuventus Plus Nr IP2010 0076/ 70 is gratefully acknowledged. ’ ADDITIONAL NOTE Originally submitted for the “Richard F. W. Bader Festschrift”, published as the November 17, 2011, issue of J. Phys. Chem. A (Vol. 115, No. 45). ’ REFERENCES (1) Gresh, N.; Cisneros, G. A.; Darden, T. A.; Piquemal, J.-P. J. Chem. Theory Comput. 2007, 3 (6), 1960–1986. (2) Engkvist, O.; Astrand, P.-O.; Karlstroem, G. Chem. Rev. 2000, 100, 4087–4108. (3) Stone, A. J. Program ORIENT, http://www-stone.ch.cam.ac.uk/ programs.html. (4) Dominiak, P. M.; Volkov, A.; Dominiak, A. P.; Jarzembska, K. N.; Coppens, P. Acta Crystallogr., Sect. D 2009, 65, 485–499. (5) Dominiak, P. M.; Volkov, A.; Li, X.; Messerschmidt, M.; Coppens, P. J. Chem. Theory Comput. 2007, 3, 232–247. (6) Volkov, A.; Messerschmidt, M.; Coppens, P. Acta Crystallogr., Sect. D 2006, 63, 160–170. (7) Volkov, A.; Li, X.; Koritsanszky, T. S.; Coppens, P. J. Phys. Chem. A 2004, 108, 4283–4300. (8) Zarychta, B.; Pichon-Pesme, V.; Guillot, B.; Lecomte, C.; Jelsch Acta Crystallogr., Sect. A 2007, 63, 108–125. (9) Dittrich, B.; Hubschle, C. B.; Luger, P.; Spackman, M. A. Acta Crystallogr., Sect. D 2006, 62, 1325–1335. (10) Ba) k, J. M.; Domagaza, S.; H€ubschle, C.; Jelsch, C.; Dittrich, B.; Dominiak, P. M. Acta Crystallogr., Sect. A 2011, 67, 141–53. (11) Ayers, P. W. J. Chem. Phys. 2000, 113, 10886–10898. (12) Hansen, N. K.; Coppens, P. Acta Crystallogr., Sect. A 1978, 34, 909–921. (13) Clementi, E.; Roetti, C. At. Data Nucl. Data Tables 1974, 14, 177–478. (14) Hirshfeld, F. L. Theor. Chim. Acta 1977, 44, 129–138. (15) Bader, R. F. W. Atoms in Molecules: A Quantum Theory; Clarendon Press: Oxford, U.K., 1990. (16) Bader, R. F. W. J. Am. Chem. Soc. 1979, 101, 1389. (17) Bader, R. F. W.; Nguyen-Dang, T. T. Adv. Quantum Chem. 1981, 14, 63. (18) Jackson, J. D. Classical Electrodynamics; John Wiley and Sons, Inc.: New York, U.S.A, 1975. (19) Edmonds, A. R. Angular Momentum in Quantum Mechanics; Princeton University Press: New Jersey, U.S.A., 1974. (20) Sokalski, W. A.; Maruszewski, K.; Hariharan, P. C.; Kaufman, J. J. Int. J. Quantum Chem. 1989, 14, 111–131. (21) Ke) dzierski, P.; Sokalski, W. A. J. Comput. Chem. 2001, 22, 1082–1097. (22) Rafat, M.; Shaik, M.; Popelier, P. J. Phys. Chem. A 2006, 110, 13578–13583. (23) Allen, F. H Acta Crystallogr., Sect. B 2002, 58, 380–388. 1546

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547

The Journal of Physical Chemistry A

ARTICLE

(24) Wilson, A. J. C.; Prince, E. International Tables for Crystallography; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1992; Vol. C. (25) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (26) Becke., A. D. J. Chem. Phys. 1993, 98, 5648–5652. (27) Lee, C.; Yang, W.; Parr, R. G Phys. Rev. B 1988, 37, 785–789. (28) Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta. 1973, 28, 213–222. (29) Koritsanszky, T. S.; Volkov, A. V.; Coppens, P. Acta Crystallogr., Sect. A 2002, 58, 464–472. (30) Volkov, A.; Macchi, P.; Farrugia, L. J.; Gatti, C.; Mallinson, P.; Richter, T.; Koritsanszky, T. XD2006, A computer program for multipole refinement, topological analysis of charge densities and evaluation of intermolecular interaction energies from experimental or theoretical structure factors, version 5.36.; 2006. http://xd.chem.buffalo.edu/ (31) Clementi, E.; Raimondi, D. L. J. Chem. Phys. 1963, 38, 2686–2692. (32) Macchi, P.; Coppens, P. Acta Crystallogr., Sect. A 2001, 57, 656–662. (33) Coppens, P. X-Ray Charge Densities and Chemical Bonding; Oxford University Press: New York, U.S.A., 1997. (34) Laidig, K. E. J. Phys. Chem. 1993, 97, 12760. (35) Popelier, P. L. A.; Bone, R. G. A. MORPHY98; UMIST: Manchester, England, 1998. (36) Jayatilaka, D.; Grimwood, D. J. Tonto: A Fortran Based ObjectOriented System for Quantum Chemistry and Crystallography. Computational Science
ICCS; Springer: New York, 2003; Vol. 4, pp 142
151. (37) Pichon-Pesme, V.; Lecomte, C.; Lachekar, H. J. Phys. Chem. 1995, 99, 6242–6250. (38) Volkov, A.; Coppens, P. J. Comput. Chem. 2004, 25, 921–934. (39) Dominiak, P. M.; Volkov, A.; Li, X.; Messerschmidt, M.; Coppens, P. J. Chem. Theory Comput. 2007, 3, 232–247.

1547

dx.doi.org/10.1021/jp204010v |J. Phys. Chem. A 2013, 117, 1535–1547