{'fraction_non_alphanumeric': 0.031073001887979863, 'fraction_numerical': 0.01099354940213971, 'mean_word_length': 4.957005623242737, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inverse Quantum Chemistry: Concepts and Strategies for Rational Compound Design

The rational design of molecules and materials is becoming more and more important. With the advent of powerful computer systems and sophisticated algorithms, quantum chemistry plays an important role in rational design. While traditional quantum chemical approaches predict the properties of a predefined molecular structure, the goal of inverse quantum chemistry is to find a structure featuring one or more desired properties. Herein, we review inverse quantum chemical approaches proposed so far and discuss their advantages as well as their weaknesses.


Introduction
The last decades have witnessed the fast-paced development of a wide range of quantum chemical methods [1]. Examples are highly accurate but also computationally demanding wave-function-based approaches like Coupled Cluster and Configuration Interaction [2] and density functional theory [3]. The latter allows for the description of large molecular systems consisting of hundreds of atoms -usually at reduced accuracy, however. Today, these computational approaches are firmly established in chemistry, and they are crucial for basically all areas of chemical research.
Underlying all these methods is the time-independent Schrödinger equation, where Ψ is the sought-for wave function characterizing the system under investigation, E the energy associated with this wave function, andĤ is the Hamiltonian operator. In nonrelativistic quantum chemistry, the Hamiltonian of an assembly of (point-like) atomic nuclei and electrons is given bŷ (in Hartree atomic units, in which the elementary charge, the mass of the electron, the reduced Planck constant, and 4πε 0 take a value of one -used throughout this work, unless otherwise stated). The indices I, J and i, j run over all nuclei and electrons, respectively.
M I denotes the rest mass of nucleus I, Z I its charge number, and r ij is the spatial distance between particles i and j. The Laplacian is defined as The first two terms in Eq. (2) are the kinetic energy operators of the nuclei and electrons, while the last three terms describe the electron-nucleus attraction, the electron-electron repulsion, and the nucleus-nucleus repulsion, respectively. The rather complicated Hamiltonian of Eq. (2) is usually simplified by invoking the well-known Born-Oppenheimer approximation, in which the nuclear coordinates are treated as (fixed) parameters. We may then define the electronic Hamiltonian and the corresponding electronic Schrödinger (eigenvalue) equation The molecular structure, defined as an assembly of atomic nuclei fixed in space, is a direct consequence of this approximation.
For a given assembly of atomic nuclei and electrons, the nonrelativistic electronic Hamiltonian is unequivocally defined. The methods mentioned in the beginning aim at an approximate solution of the electronic Schrödinger equation. In principle, all observables and molecular properties of interest can then be calculated from the wave function, which contains all information that can be known about the system [4]. Therefore, quantum chemical methods provide us with a wealth of information about a certain system, but only after this very system has been specified in terms of the nuclear framework. For the design of new functional molecules or materials one encounters the reverse situation, where a desired property is known and a molecule featuring this very property is searched for.
Given the "forward" direction from structure to property, extensive screening of structures is necessary to find a molecule with a predefined property (one example for this approach is drug design). It is thus highly desirable to develop computational approaches which would render such time-and cost-intensive trial-and-error procedures obsolete. In order to achieve this goal, the direction needs to be reversed, i.e., from property to structure.
Such approaches are called "inverse approaches".
Inverse problems do not only occur in the field of chemistry, but play an important role in diverse areas such as geophysics (where the composition of the Earth's mantle is infered from gravimetric measurements [5]), medical imaging (where, for example, attempts are made to clarify the internal structure of tissues by ultrasonic pressure waves [6]) and computer vision (where three-dimensional data from lower-dimensional images is constructed [7]). Accordingly, the theory of inverse problems constitutes a whole branch of mathematics with important contributions made in the first half of the 20th century by Russian researchers (see Ref. [8] for an extensive bibliography). Clearly, since inverse problems represent an important aspect in many fields of research, the study of inverse problems in quantum chemistry may well benefit from developments made in other fields [9].
As we will see, the greatest challenge an inverse chemical application faces is the huge size of chemical space, i.e., the space of molecules accessible by contemporary synthetic protocols. It has been estimated that this size is between 10 20 and 10 24 molecules [10].
However, an estimation of the total size of chemical space highly depends on the assumptions made; for example, the number of molecules having up to 30 carbon, nitrogen, oxygen, and sulfur atoms has been estimated to be larger than 10 60 [11], while the number of proteins that could theoretically exist is roughly 10 390 [12] (for an average size of 300 amino acid residues per proteins). This might be the reason for the fact that inverse quantum chemical approaches have only emerged during the past 20 years, although the basic idea has been present much earlier [13]. It is the purpose of this work to review all existing approaches in this field and to analyze their potential for future applications.
We will first review in section 2 the mathematical concepts developed for inverse Sturm- Finally, we also highlight developments from our research group in section 8 and then draw conclusions in section 10.

Inverse Spectral Theory
Inverse spectral theory is the branch of mathematics which studies what can be deduced about the structure of (differential) operators given some of their properties are known [14]. In quantum chemistry, this translates to the question whether it is possible to reconstruct the Hamiltonian from a predefined set of eigenvalues.
The Soviet-Armenian physicist Viktor Ambarzumian [15] investigated to which extent a differential equation is defined by its eigenvalue spectrum [16]. Basic inverse questionsalthough not explicitly stated as such -have been studied earlier in a purely mathematical context (see, for example, Ref. [17]). In fact, Lord Rayleigh briefly discussed an inverse problem already in 1877 in his work about acoustic waves [18]. Moreover, Ambarzumian was well aware of the fact that the ability of reconstructing the Hamiltonian from its eigenvalue spectrum would yield great insight into the structure and properties of matter [16]. He studied the simple case of an oscillating string, mathematically modeled by a differential equation of second order, with the boundary conditions where φ(x) describes the oscillation as a function of position x. q(x) is a continuous function, and α is a constant. He showed that the series of eigenvalues is only consistent with the differential equation, Eq. (6), if q(x) ≡ 0. It is easily seen that Eq. (6) resembles the Schrödinger equation for a particle in a one-dimensional box if we take α as proportional to the energy eigenvalue and q(x) as a position-dependent potential. Note, however, that the boundary conditions employed by Ambarzumian are different from the ones of the particle in a one-dimensional box. Ambarzumian studied socalled natural boundary conditions, where the first derivative of the function must vanish at the box boundaries, whereas for the particle in the box the function itself must be zero at the boundaries.
From the results of Ambarzumian, it follows that we can reconstruct the potential q(x)and concomitantly the full Hamiltonian -if we know the (complete) set of energy eigenvalues (see also Ref. [14]). For example, we can choose a potential q such that the resulting Hamiltonian features an eigenvalue spectrum of µ 2 π 2 , n 2 π 2 with n ∈ N and 1 < µ < 2.
Thus, we can selectively shift only the lowest eigenvalue while leaving all others unchanged.
In this case, the (unnormalized) wave function corresponding to the lowest eigenvalue µ 2 π 2 can be shown [14] to be given by where the denominator is defined as the so-called Wronskian of f and g. The function given in Eq. (9) is depicted in Fig. 10 a) for µ 2 = 3, and compared to the (regular) ground state wave function of an electron in a one-dimensional box, i.e., for q ′ (x) = 0.
[ Figure 1 about here.] The corresponding potential can easily be constructed numerically by inserting Eq. (9) into the Schrödinger equation for a particle in a one-dimensional box and solving for q. It is plotted in Fig. 10 b).
Several different methods to construct a Hamiltonian from a given energy spectrum have been developed. However, in 1946, the Swedish mathematician Göran Borg generalized Ambarzumian's work [19]. In his seminal paper, Borg showed that in the general case, one cannot determine the function q(x) from only one eigenvalue spectrum. Instead, one also needs the spectrum of a "complementary" eigenvalue problem (e.g., the same differential equation with different boundary conditions) in order to fully determine the corresponding differential operator. Thus, in the general case we cannot deduce the form of the Hamiltonian by predefining a suitable set of energy eigenvalues, as is done in the above example. It was pointed out that such an inversion is in general not possible, since different Hamiltonians (representing different molecules) can actually feature comparable properties [20,21]. Even if it would be possible to strictly invert the Schrödinger equation in the most general way, it can well be imagined that among the infinitely many different molecules, there will always be molecules that feature almost identical values for a given property. One can thus intuitively understand that for a given application, several molecules can exhibit the correct property. Therefore, different approaches have to be devised to tackle this difficult problem. In the following, we will review the most important of these approaches.

Inverse Perturbation Analysis
In molecular spectroscopy, one is often interested in determining the potential leading to a certain energy spectrum as accurately as possible. Nowadays, there exist many different approaches to accomplish this task [22][23][24][25], and there are also specialized computer programs available (see, e.g., Ref. [26]). For the sake of brevity, we shall not extensively dwell here on the details of all these approaches. However, we would like to highlight one method, which is very general in its scope of use and can also serve as motivation for the development of more advanced inverse approaches, as we shall see.
In 1975, Kosman and Hinze developed a method they called Inverse Perturbation Analysis (IPA) which allows to calculate accurate potential energy curves [27]. They started from the one-dimensional Schrödinger equation for a diatomic molecule (note, however, that the method itself is not limited to diatomics), where µ is the reduced mass of the molecule under consideration, r the internuclear distance, R(r) the nuclear wave function, and E the corresponding energy eigenvalue. U(r) is the potential energy function, which is usually not known exactly. However, many approximate potentials are known, which are often very close to the unknown exact potential like the well-known Rydberg-Klein-Rees potential. One of the central ideas of IPA is now to write the exact potential as a sum of a known approximate potential U (0) and a (small) correction term ∆U, If the approximate potential is a good approximation to the exact curve, the correction term is expected to be small, and can consequently be searched for with an iterative perturbative approach, where ∆U is the perturbation. First, the Eq. (11) is solved numerically for the wave functions R (0) and energy eigenvalues E (0) with the approximate potential U (0) . Since the true energy eigenvalues E can be obtained experimentally by spectroscopic techniques, one can calculate an energy correction which can be approximated by the first-order energy correction If one expands the potential correction into a set of predefined suitable basis functions, then the energy correction becomes which can be solved by standard methods to obtain the best set of coefficients c i . Once this set is obtained, one can use the coefficients to improve on the representation of the actual potential energy curve by employing Eq. (12). With this improved potential, one can iterate the procedure outlined above to obtain an even better potential function.
Thus, one can stepwise improve the description of the potential energy until the energy correction term ∆E is below a given threshold.
In their original paper, Kosman and Hinze successfully applied their method to HgH.
However, they found that their results heavily depend on the particular choice of basis functions f i (r). While Kosman and Hinze employed global polynomials (i.e., polynomials defined over the entire definition range of r, −∞ ≤ r ≤ ∞), Hamilton and coworkers relied on local Gaussian functions and found accurate and rapidly convergent results several years later [28]. Still, Hamilton and coworkers experienced some numerical instabilities (the results were very sensitive to changes in the parameters of the basis functions). These instabilities are a characteristic of inverse problems: they are typically ill-conditioned. In our case, the behavior of the potential energy function U(r) in regions where the wave function R(r) is zero is not determined, as can be seen by inspection of Eq. (14), and, therefore, the solution of this equation is not unique. We may recall that it can be mathematically proven that a single energy eigenvalue spectrum is in general not sufficient to determine the underlying potential uniquely. To remedy this problem of numerical instability, Wu and Zhang proposed to solve Eq. (14) by means of a singular value decomposition and found that their approach is fast, accurate and numerically stable [29].
We thus have a working method to produce accurate molecular potential energy functions, given the exact energy spectrum is known and a good approximate potential can be found as a starting point. This is certainly very interesting for the inverse point of view.
However, we are still confronted with two major problems: first, we are usually not that much interested in a molecule featuring a given energy eigenvalue spectrum, but rather in a molecule with a specific property (such as, e.g., dipole moment). These properties can, however, usually not be directly related to a specific energy spectrum. Even though there appears to be no general solution of the first problem mentioned above, we shall see in the sections to come that we might find other ways of directly relating the molecular potential to properties other than the energy eigenvalue spectrum, thereby making this approach highly interesting. We close this section with a reference to work by Politzer and Murray [30,31] who elaborate on the central role of the electrostatic potential created by atomic nuclei and electrons as the determinant of molecular properties and reactive behavior.

Model Equations
A completely different and very interesting approach is the development of (simple) model equations. This technique has been successfully employed by Marder et al. in 1991. They describe possibilities to increase the first electronic hyperpolarizability of conjugated organic molecules [32]. We illustrate the underlying design approach in the following, but note that the design of molecules with optimal optoelectronic properties has been studied by many groups [33][34][35].
If a molecule is placed in an external electric field, its charge density is polarized. When the external electric field has only a comparatively small effect on of the electronic structure of a given molecule, the polarization P can be written as a power series where E is the external electric field, and the constants α, β, and γ are the polarizability, first hyperpolarizability, and second hyperpolarizability, respectively. We note in passing that in the general case, the (hyper-)polarizabilities are no scalar quantities but instead described by tensors. However, for the following discussion, our approximation does not constitute any loss of generality. Since the electric field constitutes only a minor perturbation of the molecule, perturbation theory is a suitable method to calculate α, β, and higher order corrections (in cases where the wave length is much larger than the size of the molecule, the perturbation Hamiltonian describes the interaction of the molecular dipole moment with the external electric field). The first hyperpolarizability β can thus be obtained by correcting the wave function to second order in the electric field [36]. The resulting expression is rather complicated, involving sums over all occupied orbitals. However, it is often observed that the first hyperpolarizability is dominated by contributions from only the first few excited states. This motivates a simplification by taking only the ground state "g" and the first excited state "e" into account. The resulting two-state approximation [36], with the matrix elements defined as The influence of the three quantities X 2 ge , X gg − X ee , and 1/ (E g − E e ) 2 on the first hyperpolarizability can be studied at a simple two-orbital system as an example [36], for which we may write the wave function as a linear combination of these two orbitals, We furthermore define the energy expectation values and such that the energy difference between the two orbitals φ 1 and φ 2 is 2∆. The coupling strength between the two orbitals is then defined as With these definitions, we can write the time-independent Schrödinger equation in matrix form as For this system of homogeneous equations to have a nontrivial solution, i.e., c 1 , c 2 = 0, the determinant of the 2×2 matrix needs to be zero, which yields an equation to determine the energy E. With the two solutions and the first equality of Eq. (24), we can establish the relation For each of the two energy eigenvalues, we get a separate set of coefficients c i , such that we obtain two wave functions, and Expressions for the individual coefficients can be obtained from the normalization condition, i.e., We can thus calculate expressions for X 2 ge , X gg − X ee , and 1/ (E g − E e ) 2 as well as for β as a function of ∆. These functions are given in Fig. 10. We see that the combination of the individual expressions leads to a maximum of β at about ∆ = 0.25 when the overlap between φ 1 and φ 2 and the coupling strength t are set to 0.5. The task is now to find a molecule corresponding to these values, which is by far not trivial.
[ in such a way that they were able to come up with specific design strategies to construct the desired molecules [32].
Unfortunately, the deduction of an actual molecular (sub-)structure is not straightforward at all, but requires often a considerable amount of chemical knowledge and intuition.
Especially this last fact is a notable drawback, limiting the applicability of this method to a rather narrow subspace of chemical compound space. Therefore, this approach suffers from the same fundamental problem as the QSPR technique (see below).

Optimized Wave Functions
This severe limitation has long been recognized [37]. In 1996, Kuhn and Beratan proposed a new strategy in order to overcome this problem, which can be regarded as a first truly inverse approach in rational compound design [37]. They suggested to search for an optimal molecular structure by optimizing the corresponding mathematical objects describing it, i.e., the Hamiltonian and its wave function. The authors also point out that this inverse problem is ill-conditioned as not enough information is present in order to carry out the optimization (see also Section 3). This problem can be lifted by introducing additional constraints. For example, bonding in molecules is local in nature, leading to characteristic patterns in the Hamiltonian matrix. This can be conveniently illustrated at the example of the simple Hückel molecular orbital method [38,39], where a linear conjugated hydrocarbon chain corresponds to a tridiagonal Hamiltonian matrix due to a nearest-neighbor approximation. The overlap of basis functions located at distant atoms is in fact usually negligibly small, and thus also in more advanced electronic structure theories the Hamiltonian matrix is usually sparse. With the Hückel method, for a linear chain of three orbitals centered at x = {0, 1, 2}, we can write the Hamiltonian in matrix form as where the α i represent the orbital energies and the β i the coupling between adjacent orbitals. Since we are free in our choice of energy origin, we can arbitrarily set one of the α i to zero. Furthermore, we can also choose any suitable unit for the energy, which allows us to set one of the β i to −1. With this, we can write The task is now to find values for the remaining parameters such that a given property is optimized. As an example, Kuhn and Beratan took the electric dipole operator, the matrix elements of which can be written as if overlap contributions are neglected (c 2 + 2c 3 . Within our simple approach, µ 23 is a function of only three variables, namely α ′ 2 , α ′ 3 , and β ′ 2 . If we arbitrarily set β ′ 2 = β 1 = −1, we can plot µ 23 as a function of only α ′ 2 and α ′ 3 , as shown in Fig. 10. It is straightforward to optimize µ 23 with respect to these two parameters. For α ′ 2 ≥ −5, one finds that µ 23 is minimized for α ′ 2 = −5.00 and α ′ 3 = 0.00, which yields µ 23 = −0.97 [37]. The global minimum, however, appears to be at α ′ 2 = −4848.59 and α ′ 3 = 0.00. With these values, one finds µ 23 = −1.00.
[ Figure 3 about here.] In their pioneering study, Kuhn and Beratan proposed also to optimize not only the elements in the Hamiltonian matrix, but also the coefficients of the wave function (while the energy eigenvalues are held fixed). They applied their method to a few very simple "toy" examples as the one given above in order to demonstrate the general concept.
However, for many years after its publication, the method could not be applied to find an actual molecule exhibiting a desired property, as was hoped by the authors. The problem was not so much the sheer size of chemical compound space, but rather the fact that it is extremely difficult to construct a molecule from its Hamiltonian matrix or wave function [40]. However, in 2006, Beratan and coworkers established a general framework which can circumvent this problem (see Section 7.3).

Quantitative Structure-Activity Relationships
Quantitative Structure-Activity Relationships (QSAR) and the closely related Quantitative Structure-Property Relationships (QSPR) are central for design attempts in chemistry, biology, and materials research [41][42][43][44][45][46][47]. The origins of QSAR trace back to as early as 1863, when Cros described a relation between the toxicity of primary aliphatic alcohols and their solubility in water [48]. However, the actual mechanistic basis of contemporary QSAR approaches was laid in 1964 by Hansch and Fujita with the development of the linear Hansch equation [49], which was based on the famous work of Hammett, who proposed an equation to relate equilibrium constants of a particular reaction involving benzene derivatives to some structural parameter, which depends on the particular type of benzene substituents, and a reaction parameter, which depends on the type of reaction under study [50,51].
The basic assumption of QSPR is that every physical, chemical, and biological property depends in a systematic way on the underlying molecular structure. The functional form of this dependence is then searched for by trying to establish a model relating a given property to one or more descriptive parameters of the structure (so called descriptors). In a typical QSPR study, a training set of molecular structures and corresponding properties is taken as input. The descriptors are then usually calculated directly from the molecular structure. As descriptors, diverse quantities such as the molecular weight, the number of ring systems, the dipole moment, or the solvent-accessible molecular surface can be employed [44]. Even though one could also use descriptors which are obtained experimentally, it is advantageous to employ descriptors which can be obtained directly from the (three-dimensional) molecular structure without the need of time-consuming measurements. Once all desired descriptors have been obtained, a statistical analysis of these data is carried out in order to find a mathematical relation between the descriptors and the target property. In some cases, a simple linear function can already be used to relate one descriptor to a target property, while in other cases more complex, nonlinear, multivariate expressions might be necessary. In contemporary QSPR studies, sophisticated methods such as artificial neural networks and support vector machines are employed in order to establish a functional dependence between properties and descriptors [46]. Once such a dependence has been established, one can easily predict the described properties for arbitrary molecules (although one has to be careful when predicting a property for a molecule which is very different from any molecule of the original training set [47]).
Let us illustrate this at a very simple example. Given a set of nine trigonal-bipyramidal molybdenum complexes shown in the top right corner of Fig. 10, and would like to find an expression giving the Mo-N binding energy as a function of some descriptor(s). As the Mo-N bond length can be expected to depend on the binding energy, it appears natural to choose this bond length as a descriptor. Fig. 10 shows the binding energies as a function of the Mo-N bond length [52]. In a first approximation, this dependence can be described by a linear function. With a least-squares fit, the binding energy E can be obtained from the Mo-N bond length d as As one can see from [ Figure 4 about here.] Of course, the development of a robust QSPR for a given quantity is not trivial, and cannot be simply automatized. In particular, it is often not a priori clear which descriptors should be chosen. Furthermore, a given QSPR is usually only valid for a comparingly small subset of chemical compound space, namely for molecules which are very similar to each other in terms of their molecular structures. Therefore, these techniques are only of limited applicability for a truly inverse approach, to which the entire chemical compound space should be accessible. Moreover, even though it is easy to identify the descriptor values necessary for a given desired property once a mathematical relation between this property and its descriptors has been found, the construction of a molecule corresponding to a given set of descriptors is not straightforward at all. An approach to tackle this problem is denoted as inverse QSAR [54,55]. In this approach, one aims at constructing a library of chemical compounds similar to a given lead compound (which is known to exhibit some favorable property). In this context, chemical similarity is defined as the Euclidean distance of two molecules in the space spanned by all descriptors used. Therefore, the more similar the descriptor values are for two compounds, the more chemically similar they are expected to be. The generation of new molecules is accomplished by connecting different predefined building blocks with each other in order to obtain a chemically reasonable structure. Then, the descriptors of this new structure are computed. They can be used in an optimization algorithm such as simulated annealing in order to find structures which exhibit a predefined similarity to the lead compound. One can thus identify further compounds worth to be closely investigated without the need of screening huge molecular libraries.

Advanced Sampling Techniques
Powerful computer systems and clever optimization algorithms (e.g., simulated annealing [56], genetic algorithms [57], and Monte Carlo sampling [58]) allow us to sample a comparingly huge part of chemical space. The application of such algorithms is the most important approach in rational compound design. Therefore, we shall also review these approaches, although they are not strictly inverse, since they do not directly aim at pre-dicting a structure featuring a predefiend set of properties, but rather follow the more traditional approach where the properties of a predefined molecular structure are calculated.

Rational Design in Solid-State Chemistry
During the early 1990s, Schön and Jansen developed a sophisticated technique for synthesis planning [59,60].

Inverse Band Structure Approach
where S i is the atom type at lattice site i. A given property P (such as the band gap) is a function of the atomic configuration. The goal of the optimization is to find a configuration with one or several predefined target properties P (target) , i.e., one attempts to minimize the function where ω α is the weight assigned to the property P α . In the original approach by Franceschetti and Zunger, this function is minimized by a simulated annealing algorithm. Atomic configurations are generated by Monte Carlos moves (e.g., changing the type of a given atom).
Even though this algorithm is able to find the global minimum significantly faster than it would be possible with exhaustive enumeration, it still requires sampling of several thousand atomic configurations. It is therefore very important to have a fast method at hand to evaluate the properties of a given atomic configuration. Franceschetti and Zunger relied on a valence-force-field method to optimize the positions of the atoms of a given atomic configuration [65], and utilized a semiempirical Hamiltonian, where each atom is replaced by a pseudopotential [66], in order to evaluate the properties. Furthermore, they applied a special technique for the diagonalization of the Hamiltonian, which focuses on the few eigenvalues around the band gap [67].
Franceschetti and Zunger predicted [64] the configuration of an Al 0.25 Ga 0.75 As alloy with a maximum band gap. Recently, Zunger and coworkers studied a quaternary (In,Ga)(As,Sb) semiconductor [68]. The original approach has continuously been further developed and improved. For example, it is also possible to employ a genetic algorithm instead of simulated annealing [69,70], which generally leads to a faster convergence.
We should note here that the general method of finding the optimal species with regard to a given property in a subset of chemical space is also employed and further developed by many other groups, although they do not specifically aim at an inverse approach.
For example, Nørskov and coworkers have identified the 20 most stable four-component alloys of 32 transition, noble and simple metals in both face-centered cubic and bodycentered cubic structures [71]. For this, they searched a subspace of no less than 192'016 compounds by means of an evolutionary algorithm [71]. Given the limited space here, we refer the reader to Refs. [72][73][74][75][76][77][78] for further examples of computational materials design and to Refs. [79,80] for an example of environment-enhanced reactivity provided by inverse solvent design of a condensed phase reaction.
The approach by Franceschetti and Zunger is specifically tailored for the design of solidstate materials, and its extension to the design of isolated molecules is not necessarily straightforward. In the above procedure, the total number of atoms is defined from the outset, as a certain unit cell with a given number of lattice sites is defined. However, for an isolated molecule, this approach is not valid. One would therefore have to abandon this constraint and introduce different ones. For example, by means of simple valence rules, one could define the total number of atoms coordinating to a given atom (i.e., the number of nearest neighbors). Furthermore, in the design of solid-state materials one can apply periodic boundary conditions, and therefore in principle study a crystal of infinite size.
For an isolated molecule, however, this is not possible, and one therefore encounters the problem of unsaturated valences.

Linear Combination of Atom-Centered Potentials
One of the more daunting problems in the field of rational design is the fact that individual molecules represent discrete points in chemical space, whereas most optimization methods are developed to work with continuous functions. In 2006, Beratan, Yang and coworkers addressed this problem [81]. They started from the well-known fact that within the framework of density functional theory any molecular system is characterized by the external potential υ(r) [82], which represents the attraction between electrons and nuclei, and the total number of electrons. They then treated this external potential as a continuous variable to be optimized. Yang et al. even presented a functional of external potentials together with a corresponding variational principle [83]. One major issue with this approach is that not all external potentials correspond to a molecule, even though all molecules map unequivocally to a given υ(r). An external potential which corresponds to a chemical structure is called C-representable. Beratan, Yang and coworkers very elegantly solved this problem of C-representability by expanding the potential υ(r) into a set of atomic potentials [81], where b is the expansion coefficient of the potential term υ

(R)
A (r), which denotes the potential of atom A at position R. This approach has been called "linear combination of atomic potentials" (LCAP), and it constrains the external potential to be necessarily Crepresentable. The optimization problem then reduces to finding the best set of expansion coefficients for these atomic potentials. It is important to note that not only individual atoms can be represented by such atomic potentials but a collection of such potentials can be contracted to represent e.g., a functional group, and can be optimized as a whole. In order to correspond to an actual molecule, the coefficients b (R) A must be either 0 or 1 denoting either the absence or the presence of the atom or functional group, respectively, represented by υ A (r). However, since the expansion coefficients are treated as continuous variables during the optimization, one can find non-integer values for them. In this case, rounding to the nearest integer is necessary.
While this approach of continuous optimization is rather new in quantum chemistry, Zunger and coworker pointed out that Bendsøe and Kikuchi had already presented a similar approach in the context of materials design [84].
In their original publication [81], Beratan, Yang and coworkers investigated as a proof of concept a rather simple problem, constructed from two sites (fixed in space), for which they searched for the best potential such that the resulting hyperpolarizability is maximized. In doing so, they started from a small library of six different chemical substituents modeled by different potentials, and investigated how the optimal potential is composed of these "atomic" potentials. They finally found that H 2 S 2 has the largest hyperpolarizability of all molecules present in the small chemical space considered (six substituents at two sites give a total of 2 6 = 64 possible molecules) which was in agreement with the results found by exhaustive enumeration. With the same methodology, Keinan et al. successfully optimized the first hyperpolarizability of porphyrin-based nonlinear optical materials [85]. d'Avezac and Zunger recently published results showing that the approach by Yang and coworkers is very efficient for optimizing the atomic configurations of alloys [86].
In recent years, the group of Yang continuously developed and improved the method further. In 2007, they implemented it into an AM1 semiempirical framework [87]. In 2008, a gradient-directed Monte Carlo approach, which is able to "jump" between chemical structures, and thus to overcome barriers between local minima, was presented [88], and the approach was also implemented within a tight-binding framework [89]. Furthermore, these authors studied different optimization procedures and found that a combination of their LCAP approach with a best-first search algorithm (BFS, an algorithm searching a graph by stepwise expanding the most promising node [90]) improves the performance of the optimization [91].
The LCAP and related approaches have reached a state of maturity which allows them to be applied to rather complex problems [85]. However, the position of the potentials is fixed in space (and only optimized in the course of a standard structure optimization), which restricts the space of chemical compounds studied. This bears the danger of overlooking a good candidate, simply because it is not in the sampled space. This problem could be remedied by employing an extremely large library and allowing potentials at many different sites. However, this would also dramatically increase the computational cost of the resulting optimization problem. Nevertheless, as computer systems become more and more powerful, the LCAP will certainly be applied to increasingly complex problems.

Alchemical Potentials
Another interesting approach to convert chemical space into a continuous space has been fully determined by the total electron number and the external potential, which itself is a functional of the nuclear charge distribution. Rational compound design can thus be described as an optimization problem within the space spanned by nuclei and electrons, which leads us to the minimization problem [21] min N el ,Z(r) where O 0 is the desired target property. A continuous optimization in the space of electrons and nuclei implies fractional electronic and nuclear charges. At this point, we should note that some controversy concerning fractional electrons emerged some years ago and is still not fully resolved (see, e.g., Ref. [92]). In practice, in von Lilienfeld's approach fractional atomic and electronic charges are generated by means of a suitably parametrized effective core potential (ECP) [93]. Therefore, there are no explicit fractional electrons occurring in the computation, but they are rather implicitly dealt with in the ECP.
In order to efficiently search the space spanned by electrons and nuclei, von Lilienfeld et al. worked out a mathematical expression for the derivative of the total electronic energies with respect to the nuclear charge distribution [21]. This has subsequently become known as the nuclear chemical potential µ n (r) (in analogy to the electronic chemical potential), which, for a discrete nuclear charge distribution and perturbation theory to the first order, is found to be the electrostatic field E (1) (r), where the integration is over the entire space. The nuclear chemical potential is a function of space and is related to the proton affinity [21]. At the position of a nucleus, it measures the tendency of this nucleus to "mutate" into a different type [21]. At these specific points, the nuclear chemical potential is usually called the alchemical potential [21].
Von Lilienfeld et al. demanded a rigorous description of chemical space [94]. They found that such a description is tightly connected to grand-canonical ensemble theory and, hence, could benefit from earlier developments. It was soon recognized [94] that molecular properties such as the potential energy represent a state function within the space of Z(r) and N el , and hence, any two points in this space can be connected by an order parameter, which is usually denoted as λ. This parameter can be used to interpolate between any two molecules (similar to the thermodynamic integration) and to study the variation of a given property as one molecule is continuously changed into the other [94][95][96][97]. In such studies, von Lilienfeld and coworkers found that the path connecting any two molecules is generally not linear [96]. This makes the prediction of properties of unknown molecules based on results of known compounds very difficult. A general method of linearizing these connection paths has yet to be found, and von Lilienfeld even offers a prize (the equivalent of an ounce of gold) to the person finding a solution to this complicated problem [93].
So far, this approach has not been applied to many problems yet. As a proof of principle

Mode-and Intensity-Tracking
Both Mode-and Intensity-Tracking are rather different from all the above mentioned approaches in the sense that they do not aim for the rational design of a molecular compound or material. However, they belong to the realm of inverse methods as they turn the usual procedure in quantum chemical calculations upside-down [99]. More specifically, Mode-Tracking [100] allows for the iterative, but direct optimization of qualitatively predefined (local of collective) molecular vibrations without the need of first calculating the entire Hessian matrix (as would be the case in the standard approach). Similarly, Intensity-Tracking enables one to iteratively converge a spectrum based on the most in-tense vibrations starting from intensity carrying structural distortions [101][102][103][104]. Thus, these two methods directly target molecular vibrations and allow one to obtain these at lower computational cost than with the traditional methodologies.
The standard approach to calculate the k-th vibrational frequency of a molecule is to solve the harmonic eigenvalue problem where the eigenvector q k is called the k-th normal mode and the eigenvalue λ k is proportional to the square of the k-th vibrational angular frequency. H (m) is the mass-weighted Hessian matrix (in Cartesian coordinates), the entries of which are given by [105] H (m) where m I and r I denote the mass and spatial position of nucleus I (α ∈ {x, y, z}), respectively. In the following, we will abandon these two-component subscripts in favor of single-component subscripts ranging from 1 to 3N written in lower case. The idea of Mode-Tracking is to selectively calculate certain normal modes without first computing the entire Hessian matrix, which is the most time-consuming step of a standard vibrational analysis [100]. It is usually straightforward to create a molecular distortion as an approximation for the desired vibration. One can therefore start with a guess mode The last equality means that σ k can be calculated by computing the directional derivative of the nuclear gradient along the collective distortion b, which can be conveniently done in a semi-numerical fashion [100]. In Mode-Tracking, the target mode is iteratively improved by a Davidson-type subspace iteration method [100]. In every iteration l, a new vector b (l) is generated, and concomitantly, also a new vector σ (l) is obtained. The i vectors generated up to iteration i can be assembled into the matrices B (i) and Σ (i) , respectively, to produce the matrix i.e., it is not necessary to calculate the full Hessian matrix H (m) . We then solve the small eigenvalue problem where the matrices ρ (i) and u (i) contain the the i-th approximations to the target eigenvalue λ k and target normal mode q k , respectively. The best approximations, which we will denote as ρ is calculated. From this residuum vector, a new basis vector b (i+1) is generated, where X (i) is a preconditioner. The convergence characteristics of the Mode-Tracking algorithm strongly depend not only on the initial guess mode, but also on the particular choice of X (i) . However, experience shows that also very simple choices for the preconditioner (e.g., a unit matrix) facilitate fast convergence [106].
The Mode-Tracking approach has successfully been applied in studies of a large [(Ph 3 PAu) 6 C] 2+ cluster [107], molecular wires (i.e., long carbon chains) [108], adsorbate vibrations [109], and tetrameric methyl lactate clusters [110]. Furthermore, in 2008, Herrmann et al. integrated the Mode-Tracking methodology into a QM/MM framework [111]. Recently, Kovyrshin and Neugebauer [112,113] presented a similar tracking methodology for the selective calculation of electronic excitations in the framework of time-dependent DFT.
Intensity-Tracking generally works with the same principle as explained above for Mode-Tracking. However, the approach does not aim to refine a given user-defined vibration, but to yield those normal modes which are associated with highest intensities (for a given vibrational spectroscopy). Therefore, the most important difference to Mode-Tracking consists in the choice of the starting guess mode. This guess depends strongly on the kind of vibrational spectroscopy for which the highest-intensity modes are to be obtained, but can generally be obtained from an eigenvalue problem defined from this particular spectroscopy. We have developed Intensity-Tracking for conventional infrared spectroscopy [103], Raman and Raman optical activity spectroscopy [104], as well as resonance Raman spectroscopy [102].
Mode-and Intensity-Tracking demonstrate that subspace iteration methods are a feasible way to iteratively solve an inverse problem -here, to optimize normal modes from predefined molecular distortions.

Gradient-driven Molecule Construction
In 2012, our group proposed a new concept for the rational design of functional molecules [114], which we have elaborated into the concept of Gradient-driven Molecule Construction (GdMC) [115]. The ultimate criterion for the stability of a given molecule is the Gibbs (or Helmholtz) free energy. However, in (chemical) processes involving energy changes much larger than the thermal energy RT the entropic (and temperature) contributions may be neglected as the change in free energy is then dominated by the change in electronic energy.
This situation simplifies the design problem as it is governed by the electronic structure of the reactants in a chemical process -which in turn is largely determined by the external electrostatic potential of atomic nuclei. The electronic energy difference is calculated from the electronic energies of the reactants. The latter are determined as stationary points on the Born-Oppenheimer potential energy surface, which is the mathematical object that relates all electronic energies calculated for fixed nuclear frameworks of a given number of atomic nuclei (and electrons). Consequently, a characteristic of stable reactants and of first-order transition structures is a vanishing first derivative of the electronic energy with respect to all nuclear coordinates, i.e., a vanishing geometry gradient. Therefore, the geometry gradient is the central property for determining the stability of a molecule. This fact is utilized in GdMC as an optimization criterion for the search of molecular scaffolds that fulfill a certain chemical purpose [114]. For instance, one may ask the question what chelate ligand may stabilize a molecular fragment, such as a transition metal binding a small inert molecule, by designing the scaffold of the chelate ligand in such a way that the geometry gradient of each nucleus in the composed system vanishes.
The structure of the fragment may be chosen such that it represents different degrees of activation as necessitated by the overall chemical process of which the designed species may be one intermediate. Clearly, such a design may be carried out for each node of a reaction network and mutual stability dependencies may need to be fulfilled, too, resulting in a highly complex optimization problem. However, different optimization approaches are possible within GdMC and we have explored some of them in our study on dinitrogen binding molybdenum complexes [115]. In a first approach, the Schrock complex, which binds molecular nitrogen (see Fig. 5), was taken as an example of an existing nitrogenbinding complex that could be used as a known reference to validate the approach. In one optimization procedure, a range of model complexes was deduced from this compound in which a central Mo-N 2 fragment was taken from the full Schrock complex but where the ligand sphere was modeled by the variable substituents R 1 and R 2 . Then, the atoms or functional groups for R 1 and R 2 were searched such that the geometry gradient on all atoms of the complex was minimal, leading to a stable complex featuring the predefined bond lengths mentioned above.
[ Figure 5 about here.] This procedure is illustrated in Fig. 5. In a first optimization cycle, different atoms (oxygen, nitrogen, carbon, and phosphorus) were placed at R 1 and all valencies are saturated with hydrogen atoms. An analysis of the resulting geometry gradient reveals that oxygen is best suited for minimizing the geometry gradient. One can further reduce the gradient by varying R 2 in further optimization cycles. We find that with an ethyl group for R 2 , our model complex features a gradient which is already very close to the reference gradient of the original Schrock complex.
While we have taken only the first steps within the GdMC approach, we expect that it can be turned into a black-box inverse as well as high-throughput screening tool for the design of molecules with well-defined structural features in one of its fragment structures.

Conclusions and Outlook
The rational design of molecules and materials remains a challenging topic, despite the many improvements made in this field over the past decades. The most advanced contemporary methods all rely on sophisticated optimization algorithms and powerful computer hardware. Therefore, future advances in these areas will automatically allow one to tackle increasingly complicated problems. The possibility to carry out quantum chemical calculations in real time [116] will further enhance our capabilities. Still, even the most advanced optimization algorithm will eventually encounter its limits. Therefore, further progress in inverse quantum chemistry will heavily depend not only on the development of efficient optimizing routines, possibly tailored to the specific problem at hand, but also on possibilities of evaluating desired molecular properties as efficiently as possible.