|Home | About | Journals | Submit | Contact Us | Français|
A procedure to derive the electrostatic potential (ESP) for dynamic charge densities obtained from structure models or maximum-entropy densities is introduced. The ESP essentially is obtained by inverse Fourier transform of the dynamic structure factors of the total charge density corresponding to the independent atom model, the multipole model or maximum-entropy densities, employing dedicated software that will be part of the BayMEM software package. Our approach is also discussed with respect to the Ewald summation method. It is argued that a meaningful ESP can only be obtained if identical thermal smearing is applied to the nuclear (positive) and electronic (negative) parts of the dynamic charge densities. The method is applied to structure models of dl-serine at three different temperatures of 20, 100 and 298 K. The ESP at locations near the atomic nuclei exhibits a drastic reduction with increasing temperature, the largest difference between the ESP from the static charge density and the ESP of the dynamic charge density being at T = 20 K. These features demonstrate that zero-point vibrations are sufficient for changing the spiky nature of the ESP at the nuclei into finite values. On 0.5 e Å−3 isosurfaces of the electron densities (taken as the molecular surface relevant to intermolecular interactions), the dynamic ESP is surprisingly similar at all temperatures, while the static ESP of a single molecule has a slightly larger range and is shifted towards positive potential values.
The electrostatic potential (ESP) is important for understanding the chemical reactivity and the atomic structure of molecules and solids. A variety of properties can be derived from the ESP, for example atomic and anionic radii, electronegativities, and energies (Politzer & Murray, 2002 ).
The ESP is most easily computed for a single molecule or finite cluster of atoms, for which a well defined electron density is available. This is the case for the static electron densities obtained by molecular quantum chemical methods (Kumar et al., 2015 ). Considerable effort has been devoted to the development of methods for calculating the ESP from the static electron density of an isolated molecule, which is described by the multipole (MP) model (Stewart, 1976 ; Hansen & Coppens, 1978 ) as it can be extracted from a crystal structure (Su & Coppens, 1992 ; Ghermani et al., 1993 ; Stewart & Craven, 1993 ; Volkov et al., 2006 ). One method of analysis comprises the consideration of the ESP on a surface enveloping the molecule. In applications to small and large molecules up to proteins, the ESP has thus been used to identify electrophilic and nucleophilic sides, to characterize hydrogen bonds, and to analyse intermolecular interactions (Du et al., 2016 ; Kalaiarasi et al., 2016 ; Kirby et al., 2014 ; Malinska & Dauter, 2016 ; Niranjana Devi et al., 2017 ; Sirohiwal et al., 2017 ; Zarychta et al., 2015 ; Zhurova et al., 2016 ).
The ESP of an infinite crystal is not uniquely defined. For example, the correct ESP requires that all finite approximations to the infinite sum pertain to electrically neutral crystals, for instance the summation needs to be performed over complete unit cells. One solution to this problem is the Ewald summation method (Ewald, 1921 ), which combines converging sums in direct and reciprocal spaces. Earlier methods of evaluating the ESP in crystals involve its computation directly from the X-ray diffraction data (Bertaut, 1952 , 1978 ; Stewart, 1979 ). These methods tend to suffer from series termination effects of the Fourier sums. In other methods the thermal averaged deformation density is used in combination with the ESP of the static independent atom model (IAM), resulting in an ESP approximately valid for static densities of crystals (Spackman & Stewart, 1981 ; Spackman & Weber, 1988 ; Brown & Spackman, 1994 ; Spackman, 2007 ; Franchini et al., 2014 ).
Tanaka et al. (2006 ) have computed the ESP from an electron density obtained by the maximum entropy method (MEM) applied to X-ray diffraction data. Owing to the dynamic character of this density, the reciprocal-space sum of structure factors converges. The nuclear contribution is computed by Ewald summation involving the positions of the nuclei (Tanaka et al., 2009 ; Fujiwara et al., 2012 ). This method results in an ESP that combines a dynamic electron density with a static nuclear density.
Here we propose a method of computation of the ESP for dynamic charge densities inside a crystal. The method basically involves the inverse Fourier transform of the structure factors of the total charge density. However, instead of experimental structure factors, the structure factors of a model are employed, which then are required up to resolutions of at least Å−1 in order to reach convergence. Presently, we can compute the dynamic structure factors of the MP model and IAM as well as of MEM densities.
where 0 is the permittivity of free space. A charge density is defined in units of elementary charge per volume as the difference between proton and electron densities,
For a collection of atoms or pseudoatoms with atomic numbers and static electron densities centred at positions , the total charge density can be expressed by a sum over all atoms in the crystal,
This charge distribution leads to the electrostatic potential for static charge distributions,
A periodic structure has unit cells, each filled with atoms, such that = . In analogy to the structure factor of electrons, the structure factor of the total charge density is defined as
where is the aspherical atomic scattering factor of atom j (Coppens, 1997 ). The total charge density of a periodic structure can then be expressed as the inverse Fourier transform of its structure factors, where the latter are defined on the nodes of the reciprocal lattice and is the volume of the unit cell:
The summation extends up to an upper limit , which is chosen to be sufficiently large for convergence to have been reached. The ESP can also be expressed as an inverse Fourier transform involving the structure factors of the total charge density,
where the term is excluded from the summation.
In the limit of arbitrarily large crystals ( and ), both expressions for converge to the same values [equations (3) and (6)]. The same is true for the direct- and reciprocal-space expressions for the ESP [equations (4) and (7)]. However, convergence is too slow.
The convergence problem has been solved by Ewald (1921 ) in what has become known as the Ewald summation method. In this method the ESP of the total charge density is obtained as the sum of direct-space and reciprocal-space contributions,
Each of the two contributions is described by a rapidly converging series,
where is the error function and the single adjustable parameter η should be chosen such that both sums rapidly converge.
It is noted that the product actually is the structure factor of a dynamic charge density, whereby each atom has been assigned the same isotropic displacement parameter . This observation leads to the conjecture that equation (7) for will converge sufficiently rapidly if is replaced by the structure factor of the dynamic charge density [compare with equation (5)],
where is the Debye–Waller factor of atom j. For example, for anisotropic, harmonic atomic displacement parameters it is
The same idea has been used to compute dynamic electron densities by Fourier inversion of dynamic structure factors (Mondal et al., 2012 ). Convergence has been demonstrated for sums that include all structure factors up to Å−1, even in cases where atomic displacement parameters basically represent zero-point vibrations (Mondal et al., 2012 ). Owing to the additional factor of , a more rapid convergence is expected for the ESP [equation (11)]. Present calculations confirm this convergence behaviour.
Within the present approach the Debye–Waller factor is responsible for convergence of the reciprocal-space summation [equation (13)]. On the other hand the Debye–Waller factor represents thermal motion of the atoms. In order to arrive at the ESP of a dynamic charge density it is thus important to employ for each contributing atom the same position and the same atomic displacement parameter for its nucleus and its electron density [equation (11)]. Non-matching values for atomic displacement parameters of nuclei and electron densities lead to an apparent ESP without a clear physical meaning. The latter function may then contain artefacts like Fourier ripples. The choice of matching atomic displacement parameters is implicit in equations (11) and (13), and it appears the logical choice for dynamic ESPs to be based on structure models, like the multipole model and the IAM. This requirement poses a challenge for dynamic electron densities that do not originate from a model, but are, for example, obtained by the MEM. Here the model best matching to the MEM electron density should be used in equation (15). These aspects are discussed in §3.2.2.
Following earlier work on the dynamic electron density and the MEM, the electron density is described by its values on a grid over the unit cell (van Smaalen et al., 2003 ; Mondal et al., 2012 ). Structure factors at scattering vectors follow by discrete Fourier transform from the electron density, employing the fast Fourier transform (FFT) algorithm. The resolution in reciprocal space is inversely proportional to the mesh of the grid in direct space. Convergence of the Fourier transform and its inverse is obtained for a mesh better than ~0.04 Å, corresponding to Å−1 (Mondal et al., 2012 ).
In a first step a dynamic model electron density is produced by the software PRIOR for an MP model or an IAM, or a MEM density is obtained by the BayMEM software (van Smaalen et al., 2003 ). Of course, gridded dynamic electron densities from other sources, including theoretical densities, can be used as well.
A new software, dESP, has been written, which employs this dynamic electron density together with atomic coordinates and displacement parameters for generating the ESP of the dynamic charge density according to equation (11). The program dESP will be part of the BayMEM suite.
The computation involves the following steps:
(1) Load the gridded dynamic electron density.
(2) Apply the FFT to produce the structure factors.
(3) For each scattering vector of the grid, calculate the thermally smeared nuclear structure factor from atomic coordinates and displacement parameters.
(7) Write the ESP to a file.
Dittrich et al. (2005 ) have published high-resolution X-ray diffraction data of dl-serine as measured at temperatures of 20, 100 and 298 K. They reported MP models based on invariom refinements (Dittrich et al., 2013 ) against all three datasets. Their results demonstrated a consistent description of the aspherical atomic electron densities, independent of temperature. Mondal et al. (2012 ) have reported an MP refinement against the 20 K data. Subsequently invariom-like refinements were performed against the 100 and 298 K data, where the MP parameters were kept fixed at their 20 K values and only positional parameters and atomic displacement parameters were refined. Mondal et al. (2012 ) employed these three structure models for a study of the effect of temperature on the dynamic electron density. It was found that within regions of bond-critical points (BCPs) static as well as dynamic densities possess surprisingly similar topological properties.
Here we have reproduced the refinement strategy of Mondal et al. (2012 ), arriving at multipole models MP(20), MP(100) and MP(298) for dl-serine at 20, 100 and 298 K, respectively. These MP models involve isotropic displacement parameters for H atoms. H-atom distances C—H and N—H were fixed to values from the invariom database (Dittrich et al., 2013 ). Then for each temperature, anisotropic displacement parameters were computed for H atoms, employing the SHADE3 server (Madsen, 2006 ). These values were introduced into the structure models and kept fixed during subsequent refinements. This procedure resulted in three more structure models, denoted as AH(20), AH(100) and AH(298), respectively.
Dynamic electron densities were computed for all six structure models by the software PRIOR (Mondal et al., 2012 ) of the BayMEM suite (van Smaalen et al., 2003 ). The ESP was obtained from these dynamic electron densities together with the corresponding structure models, employing the newly written software dESP (§3.1).
The software BayMEM was employed for the computation of dynamic electron densities according to the MEM. MEM electron densities were generated for each of the three temperatures, employing the corresponding dynamic model electron densities, MP(20), AH(20), AH(100) and AH(298), as prior densities. A value of was used to define convergence of the MEM calculations. This value was the lowest value for which the MEM densities did not exhibit spurious maxima. The resulting dynamic MEM electron densities are denoted as MEMP(20), MEMAH(20), MEMAH(100) and MEMAH(298), respectively.
The ESP was computed by dESP for each of the four MEM electron densities, employing atomic coordinates and displacement parameters from the corresponding prior for the nuclear part. In the case of MP priors these values can be expected to be close to the unknown values hidden in the MEM density. In the case of an IAM prior significant deviations might be present, which then will lead to artefacts in the ESP.
The electrostatic potential is particularly interesting in the region between atoms. For this purpose the variation of the ESP is considered on an isosurface of the electron density, which envelopes entire molecules. We have chosen the isosurface at 0.5 e Å−3. This value is higher than the highest value of the electron density between molecules, which is the electron density in intermolecular hydrogen bonds, with a maximum value of 0.32 e Å−3 in serine. And it is substantially smaller than the electron densities in covalent bonds. For visual inspection, this isosurface is provided in a pseudo-three-dimensional representation, with values of the ESP indicated by colours. Quantitative measures of the ESP have been obtained as integral properties over this surface according to Politzer et al. (2001 ) (see Table 1 ). Mathematical definitions are given in Appendix A . Generation of the graphical representations and computations of the integral properties of ESPs have been performed by the program MoleCoolQt (Hübschle & Dittrich, 2011 ).
The ESP of the static charge density at K has been analysed for the multipole model AH(20), which includes anisotropic displacement parameters for the H atoms. On the molecular surface defined by e Å−3, the static ESP of a single molecule exhibits its most negative values near the three O atoms, and its most positive values near H atoms attached to nitrogen and oxygen (Fig. 1 a). These properties of the ESP indicate the preferred interaction sides of the serine molecule. As is actually realized in the crystal, all three O atoms act as acceptors of intermolecular hydrogen bonds, and these three H atoms are part of a hydrogen bond too (Mondal et al., 2012 ). These features of the molecular ESP are in agreement with similar features of ESPs of other molecules (Kalaiarasi et al., 2016 ; Niranjana Devi et al., 2017 ; Zhurova et al., 2016 ). The ESP of the static electron density of AH(20) has also been computed for a cluster of unit cells. Its values on the 0.5 e Å−3 isosurface in the central unit cell qualitatively exhibit the same features as the molecular ESP (Fig. 1 ).
The ESP of the dynamic charge density inside the crystal at K is given in Fig. 2 (a) for the multipole model AH(20) on the 0.5 e Å−3 isosurface of the dynamic electron density. The general features of the dynamic ESP are similar to the features described above for the static cluster ESP on this isosurface. In particular, the range of values is nearly identical inside the static cluster ESP and the dynamic crystal ESP (Table 2 ). The absolute values and are different for these static and dynamic ESPs. This difference can be explained by the fact that the static ESP has been computed for a finite cluster, while the dynamic ESP pertains to the ESP inside a crystal. The static ESP of a single molecule has a wider range of values than the dynamic crystal ESP on their 0.5 e Å−3 isosurfaces. Again, this is explained by the molecular versus crystal character of the ESPs. It is noted that the static ESP is represented in Fig. 1 by a red-to-purple colour code encompassing its full range. A unified red-to-purple colour coding for a range of −0.776 to 1.364 e Å−1 has been employed for all dynamic ESPs. The value −0.776 is the lowest value and 1.364 is the highest value of the ESP on the 0.5 e Å−3 electron density isosurfaces of all of these ESP maps (Figs. 2 –6 ). This unified range enables a direct visual comparison of the ESPs obtained from different dynamic electron densities.
Major differences between static and dynamic ESPs appear close to the nuclei. Static charge densities have singularities at the positions of the nuclei (Mondal et al., 2012 ). Accordingly, the static ESPs are very large near the nuclei (Table 3 ). Any thermal motion – zero-point vibrations are sufficient – leads to smearing of the density and a dramatic reduction of the ESP in the neighbourhood of the atoms. Concomitantly, the electron density is increased within the low-density region between the molecules (Mondal et al., 2012 ), leading to more negative values of the dynamic ESP than of the static ESP in its region of lowest values (Table 3 ).
The effect of thermal motion on the charge densities and ESPs is strengthened with increasing temperature. The maximum value of the ESP is strongly reduced on going from 20 to 100 K and again on going from 100 to 298 K [compare AH(100) and AH(298) in Table 3 ].
On the other hand, the dynamic ESP is nearly independent of temperature on the 0.5 e Å−3 isosurfaces of the electron densities (Fig. 2 ). This visual impression is confirmed by the integral quantities computed according to Politzer et al. (2001 ). All quantities possess nearly equal values at the three temperatures of 20, 100 and 298 K (Table 2 ). Only a small reduction can be observed of the variances of the dynamic ESP of AH(298).
We have chosen as reference the multipole model with an anisotropic description of the H atoms. An alternative choice for the multipole model employs isotropic displacement parameters for the H atoms, which then can be varied in the structure refinement. Models MP(20), MP(100) and MP(298) lead to ESPs which are close to the ESPs of the corresponding anisotropic models (Fig. 3 and Tables 2 and 3 ). One can thus conclude that multipole models with either isotropic or anisotropic displacement parameters for H atoms lead to similar descriptions of the ESP.
Concentrating on the dynamic ESPs at K the IAM and invariom model have been considered. While the dynamic ESP of the invariom model is comparable to the dynamic ESP of the refined multipole model AH(20), the ESP of the IAM exhibits major differences. In particular, its range on the 0.5 e Å−3 isosurface is about half the range of the ESP of AH(20) (Table 2 ). This result is in agreement with findings for static densities of other molecules, which indicate that the multipole model is essential for extracting the correct ESP (Malinska & Dauter, 2016 ).
MEM electron densities are dynamic electron densities. Ideally, the MEM produces an unbiased electron density map corresponding to the data. In practice, the MEM electron density will depend to some extent on the choice of prior density (Prathapa et al., 2013 ). For computation of the ESP the MEM electron density needs to be combined with a dynamic nuclear density, for which the structure model underlying the prior density provides the natural choice. Here, we have computed MEM electron densities at K for different prior densities. The resulting dynamic ESPs possess similar properties on the 0.5 e Å−3 isosurfaces of the electron densities (Table 2 and Figs. 5 and 6 a) as well as similar global minimum and maximum values (Table 3 ). These similar features indicate that the MEM electron densities conform to the diffraction data with only a weak dependence on the prior. A topological analysis of MEM electron densities has shown a clear but weak dependence on the choice of prior of the MEM electron density at BCPs of covalent bonds (Prathapa et al., 2013 ). The latter correspond to density values typically between 2 and 3 e Å−3, which is substantially higher that the present isosurface of 0.5 e Å−3. Major differences have been found for the Laplacian at BCPs, with positive values in the case of IAM priors and negative values in the case of multipole and invariom priors.
Comparing values of the ESP based on MEM electron densities with those for models shows that the maximum close to the nuclei is nearly identical for the two types of maps. This can be understood from the fact that strongly positive values of ESPs are dominated by contributions from the positive nuclei, and this contribution is identical for the MEM ESP and the model ESP, where the model is the one employed for computation of the nuclear contributions to the MEM ESP. Minimum values of the ESP are found in regions far away from the molecules. They are about twice as negative for the MEM ESP as for the ESP from model densities (Table 3 ). This is in agreement with the additional smearing of MEM densities as compared to structure models, always resulting in a higher electron density of the MEM density within the regions of lowest density. Significant differences between MEM ESP and model ESP are also found for the ESP on the 0.5 e Å−3 isosurface (Fig. 5 ). The range is about 67% larger for MEMAH(20) than for AH(20), with correspondingly larger values for variances and average ESP (Table 2 ). These differences also reflect the additional smearing of the electron density in MEM maps. Apart from this scaling, the general features of the MEM ESP are the same as for the model ESPs and the static ESPs, with the most negative values near O atoms and the most positive values near H atoms.
The temperature dependence of the MEM ESP is most pronounced near the nuclei, where with increasing temperature the reduction of the ESP exactly follows the behaviour observed for the corresponding models, as explained above (Table 3 ). Any temperature dependence of the other values is much smaller than the difference between the MEM ESP and the model ESP. In particular, the ESPs on 0.5 e Å−3 isosurfaces are nearly identical for MEMAH(20), MEMAH(100) and MEMAH(298) (Fig. 6 and Table 2 ).
We have defined the electrostatic potential (ESP) for dynamic charge densities. A method is proposed for the computation of this dynamic ESP for multipole models and electron densities derived by the MEM. In particular, it is shown that the reciprocal-space summation defining the ESP converges sufficiently fast for dynamic charge densities, because of the presence of a factor with a Gaussian dependence on the length of the scattering vector, as it is provided by the Debye–Waller factor. This method is implemented in a new module, dESP, of the BayMEM software package (van Smaalen et al., 2003 ).
The dynamic ESP has been obtained for various models and MEM densities of dl-serine at temperatures of 20, 100 and 298 K, employing three sets of high-resolution X-ray diffraction data taken from the literature (Dittrich et al., 2005 ). It is found that at all temperatures the ESPs of all static and dynamic charge densities possess similar features on the 0.5 e Å−3 isosurfaces of the electron densities: the most negative values appear near O atoms and the most positive values appear near H atoms. These features are in agreement with the presence of intermolecular hydrogen bonds in the crystal, and they are in line with similar features observed for static ESPs of other molecules (Kalaiarasi et al., 2016 ; Niranjana Devi et al., 2017 ; Zhurova et al., 2016 ). Major differences between ESPs on these isosurfaces are found between MEM densities and model densities, with a 60% larger range of values for the MEM densities.
The ESP exhibits only a weak temperature dependence on the 0.5 e Å−3 isosurfaces of the dynamic charge densities (Table 2 ). A significant reduction of the ESP with increasing temperature is found near the nuclei, where this reduction reflects the increased smearing of the positive charge density of the relevant nucleus (Table 3 ). Large differences are found between ESPs of dynamic and static charge densities, since the zero-point thermal smearing is already sufficient to remove the spike at the nucleus in the static electron density.
CBH would like to thank Claude Lecomte and Mark Spackman for fruitful hints and discussions.
The ESP is analysed through its values on an isosurface of the electron density. Quantitative measures are obtained as integral properties, integrated over this surface (Politzer et al., 2001 ). Definitions of the properties are given in Table 1 . Here we give the mathematical formulae, used to compute these quantities. The electron density and the ESP are given by their values on a grid over the unit cell. α and β indicate the number of surface pixels with positive and negative values of the ESP, respectively, while . Then,
This work was funded by Deutsche Forschungsgemeinschaft grant SM 55/24-1.