Ab-initio calculations of spectroscopic (magnetic resonance) properties in supramolecular systems
The determination of local structural properties, intra- and intermolecular conformations of molecular systems and supramolecular assemblies has always been and still is a challenge for modern physics and chemistry. Many advanced techniques are capable of contributing to this quest, some of the most prominent being X-ray and neutron scattering , electron crystallography, infrared (IR) spectroscopy and nuclear magnetic resonance (NMR) spectroscopy.
In crystalline systems, scattering experiments can provide very accurate atomic coordinates. Many systems, however, lack long-range order, which is required for these scattering techniques and limits their applicability. Complementary to this, NMR experiments are able to probe local structure without the need of long-range order. While magnetic resonance techniques cannot provide the full structure in terms of three-dimensional atomic coordinates, the sensitivity to the local chemical environment of an atom is one of the key advantages of this method. Therefore, NMR is well suited to investigate molecular and supramolecular systems and their mechanisms of structure formation.
Figure 1:
Interconnection between experimental and computational techniques for the determination of structural details on the atomic and molecular scale. The focus of this work is the symbiosis of ab-initio molecular dynamics simulations, property calculations and the correspondingexperiments, indicated by red arrows.
It has become increasingly common to supplement the experimental data with adequate numerical simulations. The interplay of experimental and computational techniques for the determination of structural properties is outlined in figure 1(right).
Supramolecular and biomimetic systems are extremely complex, both in structure and dynamics – not to mention their functionality. In order to investigate such systems with magnetic resonance methods, a great deal of knowledge regarding the relationships between structure, atomistic dynamics and spectroscopic properties is crucial. This is particularly true if magnetic parameters such as chemical shifts or spin-spin couplings are to be used as sources of information. One of the specific of our group is the QC-code of the signatures of condensed phase packing effects both in NMR experiment and in quantum-chemical calculations is essential for a real understanding of such complex systems.
We look at microscopic systems of very different types. Molecules, crystalline and amorphous solids as well as surfaces and nanotubes are covered. Many of these systems can be studied with the symbiotic combination of solid-state NMR experiments and ab-initio calculations. The direct comparison to experimentally accessible data can provide significantly more insight into physical and chemical questions than either of these methods alone.
It is well-known that for the reliable calculation of NMR parameters, in particular for chemical shifts, the chemical environment of the considered atoms is of crucial importance. Besides the classical geometrical variables of molecules – bond distances and bond angles – there are also non-covalent parameters to which the NMR chemical shifts can provide access. The most prominent example among these are hydrogen bonds, for which an example is shown in the left part of figure 2: The 1H NMR chemical shift of the hydrogen bonding proton of the left imidazole molecule varies strongly on the intermolecular N-N distance.
Figure 2:
Illustration of the dependence of NMR chemical shifts on the chemical environment: An imidazole-dimer (left) and a benzene–H2 complex (right). For both, the calculated dependence of the ¹H NMR chemical shift spectra on the geometrical parameters are shown. Note that no covalent bonds are involved; the entire chemical shift variation is due to non-covalent interactions.
Further, there is a through-space effect that deserves special attention, shown in the right of figure 2: When a given nucleus (here: the lower H atom in the H2 molecule) approaches a neighboring benzene, its NMR chemical shift changes significantly before the electronic densities of the two molecules overlap with each other. This effect, shown in the graph on the right of figure 2, is known as the “ring-current”-effect. Both effects can be calculated with good accuracy from electronic structure methods. They are very useful to obtain information about structural parameters of complex molecules which can be compared against NMR experiments.
Highly accurate electronic structure methods on the explicitely correlated level of theory are available for the calculation of structure and magnetic resonance properties of small molecular systems, but not for the condensed phase. Density functional theory (DFT) represents a compromise between accuracy and computational efficiency, and has been proven many times to yield good agreement with experimental data. Although DFT based methods are sometimes erratic for particularly difficult cases of electronic structure, they allow for a parameter-free prediction of structural and spectroscopic properties in complex systems, be it molecules, crystals, surfaces or liquids and solutions.
Related publications:
- D. Sebastiani:
Structure-Property-Relationships of Supramolecular Systems from First Principles Calculations
Habilitation Thesis, June 2006 - Y.J. Lee, B. Bingöl, T. Murakhtina, D. Sebastiani, W.H. Meyer, G. Wegner and H.W. Spiess:
High Resolution Solid State NMR Studies of Poly(vinyl phosphonic acid) Proton Conducting Polymer: Molecular Structure and Proton Dynamics
J. Phys. Chem. B 111, 9711-9721 (2007) - G. Brunklaus, A. Koch, D. Sebastiani and H.W. Spiess:
Selectivity of Guest-Host Interactions in Self-assembled Hydrogen-bonded Nanostructures Observed by NMR
Phys. Chem. Chem. Phys. 9, 4545-4551 (2007) - M. R. Hansen, S. Sekharan, R. Graf and D. Sebastiani:
Columnar Packing Motifs of Functionalized Perylene Derivatives: Local Molecular Ordering Despite Long-Range Disorder
J. , submitted (2009)
Hybrid Quantum/Classical (QM/MM) Modelling
Since the 1960s, various methods for the theoretical determination of the conformations of small and large molecules have been developed. On one hand, there are classical methods that use empirical parameters. These approaches can handle up to several hundred thousands of particles - within a limited accuracy and flexibility, since they require knowledge about the types of all covalent bonds and all atomic charges in the system. On the other hand, there are more accurate and more flexible methods based on quantum mechanics (QM-methods). Due to the higher computational cost of these approaches, they can only be applied to systems of up to 500 atoms - at least with the more and more common parallel computers.
In many cases, experimentally relevant systems consist of many more particles than can be treated with QM-methods, but the accuracy of classical methods is often not sufficient. Is there a solution to this dilemma?
The concept of hybrid quantum mechanics/mechanical modeling (QM/MM) simulations has been developed for these purposes. A QM/MM calculation splits the system into two different parts. One contains a small central region of the system, treated on the quantum level and thus yielding an accurate description of the atomic structure and dynamics. The second one, which contains the big ``rest'' of the system, is modeled classically, with pre-defined force-fields between the atoms, and may therefore contain a huge number of atoms.
This approach offers a combination between good accuracy and large system size. For example, small molecules which are solvated in a large water environment can be described with the QM/MM technique. In our group, we have started to investigate organo-ruthenium compounds (which are an alternative to traditional anticancer drugs).
In cases where the central important molecule is already very large on its own, it may be neccessary to cut covalent bonds in order to make the QM sub-system small enough. In this case, particular attention has to be paid to the distortion of the electronic structure of the quantum part of the molecule due to the bond cleavage. We have developed a particular technique to minimize this distortion by fine-tuning special effective core potentials for bond-capping.
A particular research focus in our group is the calculation of NMR chemical shifts within the QM/MM scheme. After the investigation of effects of hydrogen bonding in QM/MM, we concentrate now on the improvement of the accuracy of NMR properties with respect to the border between QM and MM. Finally, the method will be applied to systems like the above mentioned organo-ruthenium compounds. Our aim, however, is the generalization of the QM/MM parameters in a way that the method can easily be transferred to various systems.
Related publications:
- O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger and D. Sebastiani:
Variational optimization of effective atom centered potentials for molecular properties
J. Chem. Phys. 122, 014113 (2005) - D. Sebastiani and U. Rothlisberger:
Nuclear Magnetic Resonance (NMR) Chemical Shifts From Hybrid DFT QM/MM Calculations
J. Phys. Chem. B, 108, 2807-2815 (2004) - S. Komin, C. Gossens, I. Tavernelli, U. Rothlisberger and D. Sebastiani:
NMR solvent shifts of adenine in aqueous solution from hybrid QM/MM molecular dynamics simulations
J. Phys. Chem. B , 111, 5225-5232 (2007) - U. Röhrig and D. Sebastiani:
NMR Chemical Shifts of the Rhodopsin Chromophore in the Dark State and in Bathorhodopsin: a Hybrid QM/MM Molecular Dynamics Study
J. Phys. Chem. B, 112, 1267-1274 (2008) - S. Komin and D. Sebastiani:
Optimization of capping potentials for spectroscopic parameters in hybrid QM/MM calculations
J. , , submitted (2009)
Finite Temperature Effects in Condensed-Phase Spectra
Since we are an ab-initio group working within an experimental - or precisely spectroscopic - environment, we are interested in the calculation of spectroscopic properties. These are especially NMR chemical shifts and nuclear quadrupole coupling constants (NQCCs), but also IR- or UV-spectra. We are trying to reproduce and predict the experimental results - especially in the condensed phase.
Experiments are usually carried out at ambient conditions, i.e. 300 K, sometimes slightly above or below this temperature. Unfortunately, most theoretical calculations use static models - physics at zero Kelvin! Quite often these models turn out to be sufficient, but in some cases the deviations between the experiment in the computer and the real one are significant. Therefore we combine Molecular Dynamics (MD) with the calculation of spectroscopic properties and try to mimic the experimental situation. MD simulations are the most demanding simulations (especially in terms of computing time), so this procedure requires the usage of massively parallel codes like CPMD or CP2K and is still only feasible for relatively small systems - since long simulations times are necessary to allow a decent sampling of the available phase space. Especially the simulation of finite temperature effects on so called ``second-order properties'', e.g. NMR chemical shieldings, are in the need of large computer clusters and long computing times.
The figure on the right shows the computed temperature dependency of the NQCC of the hydrogen bonded deuteron in benzoic acid. The dotted line represents the experimental value at ambient conditions. Our calculations revealed a huge temperature effect that has not been expected before - and which is of great interest for both experimentalists and theorists since it shows, that a direct mapping of measured quantitites onto static properties as bond lengths is not always meaningful. Theory and especially computer simulations can bridge the gap between experiment and structural quantities. Our group plans on applying the method (MD/spectroscopy) also to larger systems and second-order properties - maybe including quantum effects on the nuclei.
Related publications:
- J. Schmidt and D. Sebastiani:
Anomalous temperature dependence of nuclear quadrupole interactions in strongly hydrogen bonded systems from First Principles
J. Chem. Phys. 123, 074501 (2005) - J. Schmidt, J. Hutter, HW. Spiess and D. Sebastiani:
Beyond Isotropic Tumbling Models: Nuclear Spin Relaxation in Liquids from First Principles
ChemPhysChem 9, 2313-2316 (2008)
Adsorption of small molecules on transition metal surfaces
This project is done in collaboration with the Multiscaling project in the Kremer group, in particular with Luigi Delle Site.
Transition metal surfaces play a crucial role for many reactions in heterogeneous catalysis. Their catalytic functionality can be affected by a variety of factors, such as the morphology of the surface, defects, or poisoning. The adsorption of small molecules from environment can significantly modify the catalytic efficiency of such metal surfaces.
One particular case in this scenario is the adsorption of water. From the view of an adsorbing water molecule, the surface has to compete thermodynamically with larger water clusters or simply the gas phase. Both phases provide a significantly larger entropic contribution to the free energy, which has to be compensated by a corresponding energy difference. Therefore, the theoretical investigation of the structural and energetic properties of the initial adsorption process on realistically modeled surfaces deserves particular attention.
The electron density plots compare the adsorbed systems with the water oligomer clusters and the isolated surfaces, the scale is given in units of e/A2. The formation of a weak bond between the surface nickel atom and the oxygen is clearly visible through the displacement in electronic density (dark green regions). Most of the electronic density is taken from the bonding nickel atom, which is strongly polarized, and its first neighbors. Similar to the case of adsorbed benzene, the polarization of the nickel atom which is bonded to the water molecule is significantly stronger on the stepped surface than on the flat one. This effect also translates into the higher adsorption energy, 0.403 eV for a water monomer on a stepped surface comparing with 0.242 eV on a flat one.
On both the flat and the stepped surfaces, the adsorption energies are about twice as large as it would be expected for a standard hydrogen bond (0.22 eV). Thus, the dimer adsorption on the metal surface can energetically compete with the solvation of the second molecule in liquid water, even though the optimized cluster on the surface is not directly comparable to the situation in liquid water due to the high dynamics of the hydrogen bond network at finite temperature. The electron density plots reveal the region of strongly increased electron density leading to a Ni-O bond, and an additional charge increase on top of the first water molecule. Furthermore, the amount of electronic charge density which is found on top of the bonding Ni atom is significantly stronger than for the water monomer (deep blue color compared to light blue one). It is interesting to note that the second water molecule does not bind directly to the surface, it is even repelled from it. The second oxygen is not accumulating any electronic density towards the metal surface, and there is a distinguishable region of decreased electron density (yellow color coding) below the hydrogen which points towards the surface. In contrast to this, the hydrogen bond between the two water molecules becomes slightly stronger than in the isolated dimer, as seen by the polarization of the H-bond accepting oxygen.
The third water increases the value of the total adsorption energy by almost 0.6eV, thus practically doubling the value of the dimer. It is interesting to note that this enhanced binding of the first water is practically not related to any of its geometric properties, but is rather due to the mere presence of secondary water molecules, which constitute a kind of a first solvation shell. These additional water molecules perturb and repel the electron density at the metal surface in the neighborhood of the initial molecule in such a way that the polarization of the bonding nickel is significantly increased.
In principle, the adsorption energies can be determined experimentally, but spectroscopic parameters are often easier to obtain. First-principles calculations of experimentally accessible spectra are very scarce, also because of the relatively high computational cost involved in realistic and accurate calculations. In this project, we want to bridge the gap between experiment and theory by providing ab initio calculations of IR peaks as a function of adsorption sites and cluster sizes, enabling for the first time a direct comparison of measurements and calculations.
We investigate the changes in IR frequencies when a second and third molecule are added to the initially attached one. Particular features are the modification of certain water modes due to hydrogen bonding and the appearance of new vibrational modes involving the nickel-oxygen bond. The computed red shift due to adsorption is in a qualitative agreement with experiment.
Although the numerical accuracy is not yet comparable to experiment, we show that our calculations allow us to characterize and distinguish between relevant adsorption configurations - e.g. whether the water is adsorbed on a flat surface or at a step defect. We show that a step, the simplest possible surface defect modeled by a (221) , influences significantly the vibrational modes and frequencies of the small absorbed clusters as well as the shift in frequencies due to the adsorption.
An interesting difference arises in Ni-H modes. The protons, which are not bonding to the surface may oscillate symmetrically or asymmetrically towards the surface. On flat Ni these two modes differ by 120 cm-1 while they are almost indistinguishably on the step, where the distance to the nearest Ni atom is larger for both protons.
Thus, the analysis of the initial steps of water adsorption on different nickel surfaces by means of their harmonic frequencies allows us to characterize vibrational modes and distinguish between relevant adsorption sites and sizes of adsorbed clusters.
Related publications:
- L. Delle Site and D. Sebastiani:
Effect of a step defect on the adsorption of benzene on a (221) surface of nickel
Phys. Rev. B 70, 115401 (2004) - D. Sebastiani and L. Delle Site:
Adsorption of water molecules on flat and stepped nickel surfaces from first principles
J. Chem. Theory Comp. 1, 78-82 (2005) - T. Murakhtina, L. Delle Site and D. Sebastiani:
Vibrational frequencies of water adsorbed on (111) and (221) nickel surfaces from first principles calculations
ChemPhysChem, 7, 1215-1219 (2006)
Car-Parrinello simulations of liquids and solutions
The solvation of ions in aqueous solution is of fundamental importance in chemical and biological processes. The structure and dynamics of the solvation structure of water molecules around H+ and Cl- ions can be defined via computer simulations.
Ab-initio molecular dynamics simulations represent a powerful tool to look into such structures of solvation. The approach can be directly benchmarked by comparison with experimental spectroscopic methods. Especially proton NMR chemical shifts are very sensitive to the character of the hydrogen bond network and accessible via fully periodic ab-initio calculations.
Here, the computed ¹H NMR chemical shift histograms during ~10ps Car-Parrinello MD simulations illustrate the instantaneous chemical shift distributions of two acidic samples and pure water: c=4,9M (left), c=2,6M (middle) and c=0 (right). The average NMR lines of the HCl solutions are shown as solid lines, together with a pure water reference (dotted) as a guide for the eye. The chemical shifts are referenced to pure water average value: δHH2O = 0 ppm.
To clarify the origin of the obtained chemical shift distributions, we decompose the set of ab initio proton shifts into contributions from different types of geometrical configurations: Eigen and Zundel complexes, the first solvation shells of the Cl- ions, and the regular water molecules. Here, upper picture corresponds c=2,6M, lower one c=4,9M.
The correlation between the structural arrangement in which a proton is found and the resulting chemical shifts is striking. While the regular water molecules are almost completely unaffected by the presence of the chlorine and hydronium ions, the protons which are part of an H3O+ or H5O2+ cations are significantly high-frequency shifted, with respect to regular water. The opposite is found for hydrogens in water molecules which solvate the chlorine anions: their NMR signal is located at low-frequency values.
Qualitatively, these trends can be explained via electron density considerations. In a hydronium ion, the water electrons are shared by three protons, so that with respect to neutral water, there is less electronic density available per hydrogen. The same -- to a lesser extent -- is the case for a Zundel complex. In contrast to this, the hydrogens in a chlorine solvation shell point into a large electron cloud. In addition, this cloud is somewhat more delocalized than the lone pairs of a water oxygen, due to the negative charge of the chlorine. This results in a stronger shielding of the proton spin, and hence a more negative shift than for a normal proton that is H-bonded to a water oxygen.
To validate our ab-initio results we compare the computed average chemical shifts with experimental spectra obtained by liquid-state ¹H NMR measurements for the same concentration values. The picture on the left shows a linear dependence of the average chemical shifts via acid concentration.
The very good agreement with experiment is a strong confirmation of dynamic H-bonding structure obtained from Car-Parrinello MD and ab initio NMR sampling under periodic boundary conditions. That means, the first principles calculations do give a realistic picture of microscopic structure of highly fluctuating liquids and solutions.
Related publications:
- T. Murakhtina, J. Heuft, J.E. Meijer and D. Sebastiani:
Ab-initio and experimental ¹H NMR signatures of solvated ions: the case of HCl(aq)
ChemPhysChem 7, 2578-2584 (2006) - S. Komin, C. Gossens, I. Tavernelli, U. Rothlisberger and D. Sebastiani:
NMR solvent shifts of adenine in aqueous solution from hybrid QM/MM molecular dynamics simulations
J. Phys. Chem. B , 111, 5225-5232 (2007)
Ab-initio path integral molecular dynamics and momentum densities
Quantum-nuclei-simulation isomorphic to P coupled classical systems
The path integral formalism represents an isomorphism between a quantum system and an equivalent classical model system. In the latter, each original quantum particle is represented as an ensemble of P ∈ IN classical particles (beads), which are connected among themselves with springs in a ring topology, and interact with their corresponding beads of the other (“quantum”) particles via their physical interaction potential. Due to this analogy, a particle in the isomorphic classical system is often referred to as a “ring polymer” or “ring necklace”, and the system of corresponding beads is called replica.
The conventional path integral formulation is based on the real space representation of the hightemperature density matrix. In neutron scattering experiments, however, the momentum space density n(k) = |Ψ(k)|2² is probed, which cannot be obtained directly from the regular density. Instead, the standard scheme has to be modified in such a way that the polymer is opened.
Path integrals were made popular by R.Feynman, implemented (in combination with classical potentials) and applied to superfluid helium by D.Ceperley and more recently used to investigate the quantum nature of protons in complex systems (in particular liquid water) within a density functional theory description by D.Marx.
With the help of the Trotter expansion, which corresponds to transform the density matrix to a higher temperature, it is possible to write the partition function as
with P ∈ IN.
The new aspect is that with a modification of the conventional path integral scheme, it is possible to express not only quantities in real space (R-space), but also momentum densities. A complete derivation of this modified path integral formalism would exceed the space available here, but it can be shown that the momentum density of a nucleus i, n(ki) = |Ψ(Ri)|² can be expressed as: 
using the off-diagonal direct-space density matrix n(R1,R´1):
While the conventional classical isomorphism corresponds to a ring polymer, this new modification describes a linear polymer, in which one real-space point (R1) is duplicated (yielding R1 and R´1), with new harmonic potential between these new coordinates.
Ferroelectric Phase Transitions via Path Integral Simulations
KDP or potassium dihydrogen phosphate shows a ferroelectric phase transition from an orthorhombic ferroelectric structure to a tetragonal paraelectric one at 122 K. The structure in either phase can be described as consisting of phosphate units interconnected by a H-bonding network. The body-centered tetragonal structure has PO4− units with equal P-O bond-lengths and the K+ ions are separated from the P center along the c-axis by exactly c/2.
Each oxygen in KDP is linked to another by a hydrogen-bond. Above the tranistion temperature experiments show a 50 % probabilty of finding the proton bonded to either oxygen. Below the transtion temperature all protons order to the same site (e.g., closer to all the lower oxygens of a PO4− unit along the c axis). This leads to a symmetry breaking and upper and lower oxygens (along the c axis) become inequivalent. Simultaneously, the K and P atoms also move off their symmetric sites and this cause a spontaneous electric polarization in the crsytal. It is unclear, however, which phenomenon - the proton off-centering or the K-P polarization - triggers the other. It's a classical hen-and-egg problem since both symmetry-breakings support and amplify each other.
Neutron Compton-scattering experiments on KDP also provide evidence for proton tunneling in the paraelectric phase. This tunneling process and its supression on cooling could be the initial trigger for the phase-transition. However, this mechanism of the phase tranistion is still under debate and the importance of tunneling stabilisation of the paraelectric phase as yet unaccessed.
On deuteration the transition temperature of KDP increases by 100 K, which is about a factor of two. This drastic increase has been attributed to a sum of two effects: increase in the mass of the tunneling species and the change in its potential energy profile. The latter effect appears as an increase in the oxygen-oxygen distance and is termed as the geometric effect of deuteration on the ferroelectric transition. These effects of deuteration are key in understanding the nature and mechanism of the phase tranistion in KDP.
Simulations can provide a great deal of microscopic insight into the mechanism of ferroelectric phase transition. For KDP it is essential to take into account the quantum behaviour of the nuclei (protons in particular) to allow for tunneling stabilization of the paraelectric phase. Moreover, accurate potential energy surfaces are required to simulate the phase transition realistically. Quantum nature of the nuclei are included in the statistics via Feynman's path-integral formulation of quantum statistical mehcanics.
This method essentially maps the problem of a quantum particle into one of a classical ring polymer with beads that interact through spring forces. As a result one can perform classical molecular dynamics (MD) or Monte Carlo simulations on the polymer and yet obtain the quantum statistical averages on the physical system. The implentation of path integral MD (PIMD) in the CPMD package allows this sampling to be performed on accurate (electronic) ground state potential energy surfaces.
While mapping to ring polymers (or closed PIMD) allows one to obtain averages of position-dependent operators (via configurational densities) to obtain the diagonal part of the density matrix one has to open up these polymers. This is useful for computing averages of momentum-dependent operators, in particular the momentum density distribution. Such an open-path or linear polymer scheme has been implemented recently in CPMD by our group and can be used to reproduce, for instance, the NCS experiments on KDP. In fact, we believe such momentum-density plots can be very helpful in studying various hydrogen-bonded systems where the presence of quantum mechanical tunneling is suspected.
Currently, we are using this scheme to study the ferroelectric phase transition in KDP in order to assess the importance of tunneling in the system. Similar simulations on DKDP are also being performed in order to distinguish the contributions of tunneling and the geometric effects to the huge isotope effect obseved on the phase transition.
Related publications:
- J.A. Morrone, V. Srinivasan, D. Sebastiani and R. Car:
The Proton Momentum Distribution in Water: An Open Path Integral Molecular Dynamics Study
J. Chem. Phys. 126, 234504 (2007)
Nucleus independent chemical shift maps in condensed phases
The NMR chemical shift is derived from the Larmor frequency of the nuclear spin of an atom, which describes the precession of the spin when the system is placed in a magnetic field. Since the electrons also react to the external field, the total magnetic field responsible for this precession is the superposition of the external one and the field induced by the electronic response. The nuclear shielding tensor is the negative proportionality factor between the electronically induced magnetic field, taken at the atomic position, and the externally applied one.
The induced field and thus the chemical shift tensor, however, are well-defined in all points of space, not only at the positions of the nuclei. This generalization of the chemical shift into a such a scalar field is called Nucleus Independent Chemical Shift (NICS). Using CPMD, the calculation of nucleus independent chemical shifts is done using a periodic plane-wave basis, i.e. in reciprocal space. All physical quantities are stored in Fourier space by means of their plane wave coefficients, which are defined as a three-dimensional grid of reciprocal space vectors. Operators which are defined in direct space (r-representation), like the position operator itself or the exchange-correlation potential, are applied by means of forth-and-back Fast Fourier transformations (FFT) of the orbitals. Other operators, such as the Coulomb propagator or the Biot-Savart law for transforming a current into a magnetic field, can be expressed and computed very efficiently in reciprocal space (G-representation). In this approach, a single Fourier transform is enough to obtain the desired NICS field in all points of space. In contrast to direct-space based methods, it is not necessary to specify particular points for which the NICS is to be computed. Instead, it is available after a regular NMR calculation on all points of the direct space mesh at essentially no additional cost.
NICS maps are illustrated on the example of carbon nanotubes (CNT): On the right, there are three (11,0) CNTs, which show electronic currents and induced magnetic fields:
On the left, the electronic current density along the tube (x-direction=horizontal in the picture) for an external magnetic field that is oriented orthogonal to the plot – and thus orthogonal to the tube axis. The calculation is done for an infinite periodic CNT, which causes an electronic current that is not “closed” at any finite position. On the top of the image, the dark blue region illustrates a steady current “to the right” (depending on the direction of the external field), while on the bottom of the tube, the current flows in the opposite (“to the left”) direction.
Interestingly, inside the CNT, there is a smaller current density of opposite orientation! This phenomenon is known from other aromatic systems such as benzene, where a small inner ring current is oriented in the opposite direction as the so-called “π -ring-current” around the benzene ring, which is due to the delocalized π-type orbitals.
On the right, there is the component of the induced field in x-direction (horizontal in the picture, orthogonal to the CNT) for the same orientation of the external magnetic field (which is now also horizontal in the plot – still orthogonal to the tube axis). The orientation of the ring currents shown in the above image would now be orthogonal to the image (along the tube axis), on the upper and lower border of the CNT. The effect of these strong ring currents is evident: the induced field looks similar to that of a dipole. Interestingly, the field inside the nanotube is almost homogeneous.
Here, we have the trace of the induced field tensor (the isotropic NICS map) which is the average of the previous plot for the three possible orientations of the external magnetic field. In this plot, the symmetry of the CNT becomes visible: Since there is no preferred direction outside the tube, there is almost no sign of the strong ring currents. However, inside the CNT, their influence is very pronounced! The remarkably homogeneous NICS shift exceeds everything known from isolated molecules.
These illustrations show some of the insight into the nature of the delocalized electrons in carbon nanotubes that can be obtained via the analysis of the nucleus independent chemical shift maps in periodic systems. Several other molecular systems have been studied with this technique as well, for example hydrogen-bonded calix-hydroquinone nanotubes and graphene sheets. In the future, we want to apply this technique also to other condensed-phase systems, such as amino acid crystals.
Related publications:
- B. Kirchner and D. Sebastiani:
Visualizing Degrees of Aromaticity for different Barbaralane Systems
J. Phys. Chem. A 108, 11728 (2004) - D. Sebastiani:
Current Densities and Nucleus Independent Chemical Shift Maps (NICS) from Reciprocal Space Density Functional Perturbation Theory Calculations
ChemPhysChem 7, 164-175 (2006) - M. Takase, V. Enkelmann, D. Sebastiani, M. Baumgarten and K. Müllen:
Annularly-fused Hexapyrrolohexaazacoronenes: A Multiple Interior Nitrogen-containing pi-System with Stable Oxidation States
Angew. Chem. Intl. Ed.45, 5524-5527 (2007) - D. Sebastiani and K. N. Kudin:
Response Properties of Carbon Nanotubes in Magnetic Fields
ACS Nano 2, 661-668 (2008)
Methodic developments in density functional theory calculations
Density functional theory is a powerful tool for the description of molecular and supramolecular systems of up to a few hundred atoms. Isolated molecules, periodic (crystalline) and amorphous solids, liquids and solutions can be treated efficiently with this technique, in particular in combination with molecular dynamics techniques (e.g. within the Car-Parrinello-scheme). However, the computational effort in describing such a supramolecular systems depends unfavorably on the system size. Doubling the system will typically result in an eightfold higher computational cost. The first development presented here goes in the direction of linear scaling methods, which aim at reducing this dependency.
The idea of this work is a new way of calculating approximate but accurate total energies within the framework of density functional theory. Our technique is based on an expansion of the energy functional to second order and does not require self-consistent iterations of the total density. The functional can be minimized by using the same techniques as developed for variational density functional perturbation theory. The method is ideally suited to systems composed of weakly interacting fragments, but it can also be applied to semiconductors and insulators. As a test case, the method is applied to a water dimer.
The approach is best explained using a system composed of two molecules – the water dimer, for example. We start from the Kohn-Sham orbitals of the full system, which are written as the sum of the orbitals of the isolated molecule(s) and a correction: ψi(r) = ψi(0) (r) + ψi(1) (r)
With the help of this, the energy functional may be rewritten as:
The exact meaning of the various elements of this expression can be taken from the first publication cited below. The resulting total energy for the combined system – the water dimer – is shown in the following graphics:
The triangles show the energy obtained by the simple combination of the orbitals of the unperturbed isolated water molecules, the squares show the corresponding fully optimized DFT energies, and the circles (which are actually almost identical to the squares) show the density functional perturbation theory results obtained with the technique presented above. The agreement between the fully optimized self-consistent calculation and the approximate one is excellent. This work could be the starting point for a non-self-consistent DFT-based method for electronic structure calculations, which would provide significant computational advantages compared to well-established standard techniques.
The structure of this theory, combined with the use of localized orbitals, leads itself very naturally to application to QM MM calculations. Furthermore it can also be used to construct effective polarizable potentials. Conceptually, our approach allows a more intuitive description of the interaction between molecules than a usual wave function optimization. In fact, it is possible with our technique to distinguish the distortions of the wave function of an isolated molecule caused by the interaction with others. In a fully self-consistent calculation this insight is much harder to obtain.
Another recent development in the field of density functional theory was an approximate semi-empirical way of incorporating van-der-Waals forces. We add an effective atom-centered nonlocal term to the exchange-correlation potential in order to cure the lack of London dispersion forces in standard density functional theory. The calibration of this long-range correction is performed using density functional perturbation theory and an arbitrary reference. Without any prior assignment of types and structures of molecular fragments, our corrected generalized gradient approximation density functional theory calculations yield correct equilibrium geometries and dissociation energies of argon-argon, benzene-benzene, graphite-graphite, and argon-benzene complexes. The mathematical details of this work can be found in the second article cited below.
We have investigated the results of our approach for the benzene dimer and for infinite graphene sheets. The benzene dimer served as the calibration system, because for this system, many highly accurate quantum chemical calculations exist. The question is now: How good does the approach reproduce energies for the much more complex graphite sheets? The answer is shown on the right:
These graphs show the energies of the benzene dimer in the parallel (left) and the T-shaped geometry (middle), as well as for the graphene sheets (on the right). In the inset of the left graph the additional d-channel projector is depicted in arbitrary units as a function of r, the minimum is 3.3A. Potential energy curves of the total energy of interaction are plotted as a function of the distance z. Circles correspond to normal ECPs, diamonds to the optimized ECP (OECP), and squares to MP2 data. The experimental value for the interlayer distance and energy of interaction of two graphite sheets is marked by a cross. The left hand panel shows the curve for calibration. The other two graphs represent results for the T-shaped configuration of benzene (middle panel) and for the slab of graphite (right hand panel). In the inset of the right hand panel the averaged interlayer distance z is presented during first-principles BOMD using the normal ECP (dotted line) and the OECPs (continuous line).
Related publications:
- D. Benoit, D. Sebastiani, M. Parrinello:
Accurate Total Energies without Self-Consistency
Phys. Rev. Lett. 87, 226401 (2001) - O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger and D. Sebastiani:
Optimization of Effective Atom Centered potentials for London Dispersion Forces in Density Functional Theory
Phys. Rev. Lett. 93, 153004 (2004) - O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger and D. Sebastiani:
Variational optimization of effective atom centered potentials for molecular properties
J. Chem. Phys. 122, 014113 (2005) - O. A. von Lilienfeld, I. Tavernelli, U. Rothlisberger and D. Sebastiani:
Performance of optimized atom centered potentials for weakly bonded systems using density functional theory
Phys. Rev. B 71, 195119 (2005)