We optimize the internuclear geometry and electronic structure of a model chiral system to achieve a maximal photoelectron circular dichroism (PECD) in its one-photon ionization by circularly polarized light. The electronic structure calculations are performed by the single center method, while the optimization is done using quantum alchemy employing a Taylor series expansion. Thereby, the effect of bond lengths and uncompensated charge distributions on the chiral response of the model is investigated theoretically in some detail. It is demonstrated that manipulating a chiral asymmetry of the ionic potential may enhance the dichroic parameter (i.e., the PECD) of the randomly oriented model system well beyond β1 = 25%. Furthermore, we demonstrate that quantum alchemy is applicable to PECD despite the unusually strong coupling of spatial and electronic degrees of freedom and discuss the relative impact of the individual degrees of freedom in this model system. We define the necessary conditions for the computational design of PECD for real (non-model) chiral molecules using our approach.

Photoelectron circular dichroism (PECD1–3) is a fascinating chiroptical effect of a forward–backward asymmetry in the photoemission from randomly oriented chiral molecules with respect to the propagation direction of the ionizing circularly polarized light. Predicted by Ritchie4 about 50 years ago theoretically, it was first verified in one-photon5,6 and then in multiphoton7,8 ionization experiments with chiral molecules. At the present time, it is considered to be universal with respect to the photoionization regime.9,10 Because PECD is purely an electric-dipole effect, its strength typically reaches a few percent, which is much larger than the conventional circular dichroism (CD11). As a consequence, PECD can efficiently be applied to even a sub-percent enantiomeric excess determination12 in the gas phase.1–3 Since larger PECD strengths enable higher sensitivity in the enantiomeric excess determination, several attempts have already been made to maximize this chiral asymmetry. One of the discussed routes is to optimize the properties of the ionizing light using coherent control schemes,13,14 which may yield a PECD well above 30%. Another very different route proposes fixing molecular orientation partly15,16 or even fully17,18 in space. This strategy of optimizing the properties of targets as a whole may enhance the chiral asymmetry well beyond 50%.

At present, it is understood that the magnitude of the PECD is governed by the ability of the outgoing photoelectron to probe a chiral asymmetry of the molecular potential through multiple scattering effects. Therefore, a potential route to enhance the PECD itself would be an optimization of the intrinsic properties of chiral targets to generate the strongest asymmetry of its chiral potential. For this purpose, one needs to manipulate a charge density distribution in the ion, which generates an electrostatic potential for the outgoing photoelectron. This can be achieved by changing its nuclear charges, its internuclear geometry, or the distributions of the bound electrons when, e.g., constituent atoms get exchanged or nuclear charges get merged. For instance, one may think of transforming the three terminal CH3 groups in the small chiral molecule fenchone independently into NH2, OH, or F, which already yields 43 different iso-electronic chiral modifications of this molecule from three sites alone. The fundamental question immediately arises: Which of the modifications has the strongest PECD? Unfortunately, the present knowledge of the PECD phenomenon does not allow us to answer this question without enumerating all compounds computationally.

Before performing optimization of large and real chiral molecules with plenty of degrees of freedom, it is extremely important to understand the individual impacts of different inherent properties (nuclear charges, bond lengths, electron densities, etc.) on this chiroptical effect in a systematic way. In the present work, we make the first step toward this aim. In particular, we use a minimal model chiral system with a few tuneable properties and perform theoretical optimization of the chiral asymmetry of its electronic potential to achieve the largest PECD effect. Even for this relatively small theoretical problem, performing ab initio electronic structure calculations of the PECD in the whole space of free parameters is not feasible. Therefore, we employ quantum alchemy, which relies on a Taylor series expansion19–21 and can estimate the properties of iso-electronic molecules for millions of molecules at once without enumeration. The paper is organized as follows: The present model and theoretical methods are described in Sec. II. The results of the optimization are reported and discussed in Sec. III. We conclude in Sec. IV with a brief summary and outlook.

In this section, we first discuss the presently used model chiral system (Sec. II A) and then outline theoretical methods for the electronic structure and PECD calculations (Sec. II B) and optimization (Sec. II C).

The presently used model chiral system is very similar to those introduced in our previous works.22,23 It is built of model atoms, each represented by a positive point charge Q, which simulates a net nuclear charge being partly compensated by the inner electrons bound to this atom. In addition, we introduce an electronic cloud with a negative charge q and a spherically symmetric distribution er, which is supposed to represent the outer electrons shared in the molecule. In order to be chiral, there should be at least four non-coplanar and non-equivalent atoms. Even for four atoms, there are nine independent degrees for optimization of the internuclear geometry. In the present model study, we restrict the number of geometrical optimization parameters to three. For this purpose, we place one of the atoms at the coordinate origin (Q0, q0) and three on the respective coordinate axes (Qi, qi), and assume the distances to the origin along the axis (Ri) to be free parameters. Thereby, for given distances Ri, the four atoms are always non-coplanar in a maximal way. The initial values for the geometry parameters are Rx = Ry = Rz = 3.0 a.u., which corresponds to a typical molecular bond length. Figure 1 depicts a sketch of the present chiral system and discusses quantities.

FIG. 1.

Model of a four-atomic chiral molecule used in the present work. It consists of four point charges (red dots): Q0 is located at the coordinate origin, and Qx, Qy, and Qz are on the respective axes at the distances Rx, Ry, and Rz. Each nuclear charge Q is surrounded by an electron density (blue clouds) with the total negative charge q and a spherically symmetric distribution er. For transparency, the positions and sizes of the constituent parts are depicted without keeping proportions. All parameters are given in a.u. The electron charges q0, qx, qy, and qz are kept fixed to the values given in the figure. The nuclear charges Qx, Qy, and Qz and their positions Rx, Ry, and Rz are variation parameters. Their initial values are Qx = Qy = Qz = +2.5, and Rx = Ry = Rz = 3.0. The sum of all nuclear charges is always kept fixed to iQi = +11, such that the total charge of the ion is +1.

FIG. 1.

Model of a four-atomic chiral molecule used in the present work. It consists of four point charges (red dots): Q0 is located at the coordinate origin, and Qx, Qy, and Qz are on the respective axes at the distances Rx, Ry, and Rz. Each nuclear charge Q is surrounded by an electron density (blue clouds) with the total negative charge q and a spherically symmetric distribution er. For transparency, the positions and sizes of the constituent parts are depicted without keeping proportions. All parameters are given in a.u. The electron charges q0, qx, qy, and qz are kept fixed to the values given in the figure. The nuclear charges Qx, Qy, and Qz and their positions Rx, Ry, and Rz are variation parameters. Their initial values are Qx = Qy = Qz = +2.5, and Rx = Ry = Rz = 3.0. The sum of all nuclear charges is always kept fixed to iQi = +11, such that the total charge of the ion is +1.

Close modal

In addition to geometry properties, we optimize the charge density of the model system as well. For simplicity, we fix all charges in the electronic clouds and vary only the nuclear charges. In particular, the electronic charges are fixed in a very asymmetric way to be qx = −1, qy = −2, qz = −3, and q0 = −4 (see Fig. 1). Since the model molecular ion contains ten shared electrons, iqi = −10, the total uncompensated nuclear charge must be set to iQi = +11, thereby modeling a singly charged ion with an outer electron being ionized. We assume the three nuclear charges Qx, Qy, and Qz to be free parameters and adjust the central charge Q0 accordingly. The initial values of the nuclear charges are Qx = Qy = Qz = +2.5, and Q0 = +3.5. The one-particle potential, generated by this model ion for the photoelectron, is used to solve the stationary Schrödinger equation. The first five eigenstates of the Hamiltonian are assumed to be populated by ten common bound electrons, while the sixth eigenstate represents a singly occupied HOMO of this open-shell model system, which is ionized by circularly polarized light. This eigenstate is used as the initial state of the photoelectron for the time-dependent Schrödinger equation (see Sec. II B).

The angular emission distribution of photoelectrons, emitted from randomly oriented molecules ionized by circularly polarized light, is given by the following well-known parameterization:4,24
(1)
Here, “±” stands for the positive and negative helicity of the ionizing light, PL(cos θ) are the Legendre polynomials, and the emission angle θ is defined with respect to the light propagation direction. The explicit analytic expressions for the total cross section σ, anisotropy parameter β2, and dichroic parameter β1, which describes the forward–backward asymmetry in the photoelectron emission (i.e., the PECD), can be found in previous works.25,26 Their calculation requires molecular-frame photoionization transition amplitudes Aɛℓmk for the emission of partial photoelectron continuum waves with kinetic energy ɛ, angular momentum quantum numbers ℓm, and all light polarization k = 0, ±1. One should note that the emission probability (1) must be non-negative at all emission angles θ. This yields independently for the anisotropy and dichroic parameters the following confinement intervals: −1 ≤ β2 ≤ 2 and |β1| ≤ 1, respectively. Nevertheless, when both parameters simultaneously approach their validity-interval borders, non-physical solutions with a negative emission probability at some emission angles may appear. Such situations, therefore, need to be kept in mind during the optimization (see below).

In the present study, we investigate the non-resonant energy-dependence of the anisotropy and dichroic parameters. We, thus, neglect autoionizing electronic resonances in the calculations. The respective transition amplitudes were computed by the time-dependent realization of the single center (TDSC) method,22 which was successfully applied to study PECD in the one- and multiphoton ionization of chiral molecules.27–30 The working equations of the method and details on the calculations can be found in those works. In the present case, we employed single center expansion with , |m| ≤ 15, used weak exponential laser pulses g(t)=e(tt0)2/τ2 with t0 = 3 fs and τ = 1 fs, and propagated the initial electronic wave packet for the linear and both circular polarizations during 6 fs in a radial box of 250 a.u. The pulse intensity of 1010 W/cm2 was chosen to suppress the above threshold ionization by two-photon absorption. For the initial parameters of the model chiral system discussed in Sec. II A, the binding energy of its sixth electronic eigenstate is equal to 19.94 eV. The photon energies were adjusted accordingly to release photoelectrons with kinetic energies of 4, 5, 6, 7, and 8 eV. The anisotropy and dichroic parameters, computed for this initial reference set (the first line in Table I), are depicted in Fig. 2 by black stars. As one can see in the chosen energy interval, this system exhibits a sizable PECD characterized by a dichroic parameter between about +2% and +6%.

FIG. 2.

Dichroic (a) and anisotropy (b) parameters, computed (curves with symbols) and predicted (curves, see legend) for the model chiral system with given values of the variational parameters (see Table I for data sets). All data are evaluated on integer energy values only; grids guide the eye.

FIG. 2.

Dichroic (a) and anisotropy (b) parameters, computed (curves with symbols) and predicted (curves, see legend) for the model chiral system with given values of the variational parameters (see Table I for data sets). All data are evaluated on integer energy values only; grids guide the eye.

Close modal
Optimization of the model chiral system with respect to a maximal PECD (i.e., maximal |β1| values) requires variation in the six-dimensional space (Rx, Ry, Rz, Qx, Qy, Qz), and the photoelectron kinetic energy can also be considered as an additional degree of freedom. Therefore, performing ab initio calculations of the PECD for Np values of each parameter p requires pNp calculations for each kinetic energy, which is a tremendous computational task even if Np are small numbers, especially if we want to potentially address real molecules. In order to perform an efficient screening of the variational space, we here employ a multidimensional Taylor series expansion following the quantum alchemy approach,19–21 where we approximate βi from partial derivatives at a center a in all degrees of freedom x up to and including order k using the multi-index α,
(2)
By employing Eq. (2), we assume a sufficiently analytic behavior of the anisotropy and dichroic parameters as functions of all variational parameters, including their variations with the energy across, e.g., shape resonances.23 Truncating the expansion after the first few orders drastically reduces the scaling of the computational efforts to be, e.g., quadratic in Np if all mixed partial derivatives of the second order are included, and only a few Np values are required to initiate the search.

Many libraries obtain finite difference derivatives from a fixed stencil, i.e., the set of all function evaluations entering a numerical approximation. We implemented the multidimensional multivariate arbitrary order Taylor expansion from arbitrary evenly spaced finite difference stencils in the Python package nablachem.MultiTaylor31 using findiff.32 This library also allows us to select the mixed partial derivatives to be included selectively and performs L-BFGS-B optimization over the resulting multivariate polynomial. We always choose the maximum number of available points for stencils of lower than maximal order; e.g., first order derivatives are evaluated using all available and applicable data points rather than only two. This improves the accuracy of lower derivative orders. All monomials in the Taylor expansion for which there are not sufficiently many coefficients available are assumed to be zero. In this work, this concerns terms higher than the respective order or of mixed partial derivatives including more than two different variables. Starting from the center of the series expansion, the bounded limited memory optimizer L-BFGS-B with a budget of 1000 steps was used to find extremal points within constraints on each dimension in the six- or seven-dimensional search space of β1(x), while β2 was just evaluated at the resulting point. The bounds represent the unknown trust radius of the truncated Taylor expansion, since at large changes in charges Qi or atom position Ri, instability due to the highest order polynomial term is unavoidable. While bounds for the values of βi are known (|β1| ≤ 1 and −1 ≤ β2 ≤ 2), their limits for Ri, Qi → 0 or ∞ are unknown. Because of that, Padé approximants, which often exhibit advantageous convergence properties,20 could not be used.

In order to initiate the screening procedure, ab initio calculations were performed first for the following sets of parameters: Starting from the initial reference set Qx = Qy = Qz = +2.5 and Rx = Ry = Rz = 3.0, each of the six parameters was varied independently in steps of ±0.05 and ±0.10, keeping the other five parameters unchanged. The respectively modified system was used to construct an ionic potential for the photoelectron, and, subsequently, anisotropy and dichroic parameters β1 and β2 were computed by the TDSC method for the photoelectron kinetic energies 4–8 eV. This yields 4 × 6 = 24 calculations for each of the five considered energies. In addition, in order to be able to evaluate some multivariate partial derivatives, two (out of six) parameters were simultaneously varied in steps of ±0.05 for all combinations. This yields 4×652=60 additional calculations for each energy. The results of these calculations for 85 sets of the variational parameters (including the initial reference set) are provided in the supplementary material.

These ab initio data were used to calculate the derivatives for the Taylor series in expansion (2), which was employed to screen the search space for free parameters at which the value of dichroic parameter β1 reaches its desired value. For reliability, the second order Taylor expansion was used for those predictions, even though terms of the third and fourth order expansions were accessible if they did not include partial derivatives in more than two or one variable, respectively. Even using the full information provided by the calculations, for certain relatively large variations of the bond lengths (Ri) and charges (Qi), the Taylor expansions become inaccurate, and thereby the extrapolated values of the dichroic and anisotropy parameters exceed their validity intervals. Nevertheless, before diverging for larger parameter variations, the algorithm discovered (for each considered energy between 4 and 8 eV) a local maximum of the dichroic parameter β1 that fell in the validity interval |β1| ≤ 1, and the corresponding anisotropy parameter fulfilled the condition −1 ≤ β2 ≤ 2. Those predicted values of β1 and β2 parameters were additionally checked to yield a non-negative emission probability (1) at all emission angles θ. Moreover, considering the photoelectron energy as an additional degree of freedom, we were able to extrapolate the predictions outside the chosen energy interval of 4–8 eV. The explored values of β1 and β2 were verified by the ab initio calculations using predicted values of the variational parameters (Ri) and (Qi). The respective results are discussed in Sec. III A, while a detailed analysis of the present screening procedure is given in Secs. III B and III C.

We start the optimization in the chosen photoelectron energy interval of 4–8 eV and constrain our search space to maximum changes in ΔR = ΔQ = 0.1. Here, each energy is considered independently, i.e., no partial derivatives between the photoelectron kinetic energy and other parameters have been included, which allows us to evaluate the smoothness of the predictions.

TABLE I.

Reference system (Ref) and optimization results (sets 1–3) together with their estimated (β1e) and calculated (β1c) dichroic parameters. Set 1 predicts a maximal β1e value in the energy interval 4–8 eV, set 2 predicts a maximal β1e value in the extended energy interval of 2–10 eV, and set 3 predicts a minimal β1e value in the energy interval 2–10 eV.

RxRyRzQxQyQzEβ1eβ1c
Label(a.u.)(a.u.)(a.u.)(a.u.)(a.u.)(a.u.)(eV)(%)(%)
Ref 2.5 2.5 2.5 +5.6 +5.6 
Set 1 2.9 3.1 3.1 2.6 2.6 2.6 +24.1 +16.1 
Set 2 2.7 3.074 2.95 2.6 2.6 2.6 +32.7 +25.8 
Set 3 3.1 2.9 2.9 2.4 2.4 2.6 −5.1 −14.7 
RxRyRzQxQyQzEβ1eβ1c
Label(a.u.)(a.u.)(a.u.)(a.u.)(a.u.)(a.u.)(eV)(%)(%)
Ref 2.5 2.5 2.5 +5.6 +5.6 
Set 1 2.9 3.1 3.1 2.6 2.6 2.6 +24.1 +16.1 
Set 2 2.7 3.074 2.95 2.6 2.6 2.6 +32.7 +25.8 
Set 3 3.1 2.9 2.9 2.4 2.4 2.6 −5.1 −14.7 

The values of the dichroic and anisotropy parameters predicted for the optimization result set 1 from Table I are depicted in Fig. 2 by red dashed curves. The top panel of this figure predicts a maximal value of β1e(4eV)=24.1%. However, the ab initio calculations yield a somewhat smaller value of the dichroic parameter for that configuration, and this energy β1c(4eV)=16.1% (red triangles). As one can see, the predicted values of β1 are overestimated and those for β2 are underestimated (cf. the red dashed curve and red triangles in Fig. 2), while the overall trend is well-reproduced, as exemplified by the correct identification of a turning point in β2 at about 5 eV. Following the computed trend toward smaller kinetic energies, we find a maximal value of about β1c(1.94eV)=25.3% at the kinetic energy of 1.94 eV. Here and below, we list spline-interpolated extreme values of the computed dichroic parameter at non-integer energy values. In comparison, the largest value for β1 for all ab initio data points in the finite difference stencil is 9%. Since there is no formal convergence guarantee for the expansion of many properties, including βi, this numerical evidence is important as it indicates that both the dichroic and anisotropy parameters can be optimized via perturbative treatment and are sufficiently smooth and differentiable to allow for this gradient-based optimization. Moreover, we find that the predictions in both parameters βi yield clear trends and smooth curves along the photoelectron energy, even though each energy value of 4–8 eV has been predicted independently. This makes numerical effects an unlikely cause for the overestimation of the magnitude of βi but rather points toward relevant higher order terms.

In the second step, we now consider photoelectron energy as a free parameter and search for optimal β1 values within the extended interval of 2–10 eV. To this end, we have to include mixed partial derivatives between the photoelectron energy and all other parameters, i.e., we only build one model for the whole energetic, spatial, and electronic domain covered by the stencil.

From the results of set 1, one might expect the largest prediction value for energies below 4 eV. In addition, the optimization on the Taylor expansion yields a new point, set 2 (see Table I), which is predicted to have the largest dichroic parameter value at the boundary of the search interval, at 2 eV. For this system, the perturbation expects an almost linear dependency of the dichroic parameter on the photoelectron energy; see the blue dotted curves in Fig. 2(a). For the domain of 4–8 eV, we can compare this to the individual results, where five independent models predict the same configuration in set 2 without the inclusion of any mixed partial derivatives with the energy. Doing so would yield the same linear trend, although with a steeper slope for β1 and even the wrong sign of the curvature for β2. Since either method overestimates the parameter magnitude, this points toward those mixed partial derivatives containing valuable information about the interaction between the energy parameter and the spatial and electronic degrees of freedom. This will be important to consider when optimizing βi for real molecules.

Extending the configuration of set 2 (see Table I) toward even lower photoelectron energies does not reproduce the local maximum in β1 from the ab initio calculations around 2 eV, so second order mixed derivatives in the photoelectron energies are not sufficient. Note that while fourth-order derivatives in energy alone are available (since the data set contains the reference system for five different energy values), this would only apply to a frozen spatial and electronic configuration. Once those dimensions are allowed to change, their interplay can only be captured up to second derivatives in most cases, given our particular stencil. The optimal predicted value of the dichroic parameter is β1e(2eV)=32.7% [see Fig. 2(a) and Table I], while the ab initio calculations yield a somewhat smaller value of β1c(2eV)=25.8%, and the computed dichroic parameter reaches its maximal value of β1c(1.56eV)=26.3% at the lower energy of 1.56 eV. As shown in Fig. 2(b), predictions of β2 are of comparable accuracy for both set 1 and set 2. Given the perturbative nature of the model, some error from higher order terms is expected, but trends are well-preserved in the optimization procedure: just as predicted, the value of the dichroic parameter computed for the parameter set 1 at 4 eV is indeed slightly larger than that computed for set 2, and, simultaneously, its value computed for the parameter set 2 at 2 eV is slightly larger than that computed for set 1 [see Fig. 2(a)].

Finally, we search for configurations where the forward–backward asymmetry in the photoemission changes its sign, i.e., where β1 is negative. The optimizer suggests set 3 (see Table I) as the point of minimal β1, where the predicted dichroic parameter reaches β1e(2eV)=5.1%. Ab initio calculations confirm this negative value and yield even stronger negative asymmetry of β1c(1.7eV)=14.7%. While this is in line with the optimization goal, the prediction accuracy for both parameters in this task is worse than in the previous tasks. Potentially, a smaller displacement in the finite difference stencil may improve the derivative calculation, as the noise of the ab initio calculations upon small perturbations during the validation of the prediction was found to be negligible.

As discussed in Sec. III A, the models built from 85 single point calculations can be used to optimize the PECD parameters βi. This is computationally feasible now because the models are much faster to evaluate than one ab initio calculation. We can now use these models to draw statistical conclusions about which terms typically contribute strongly to the values of βi: All terms that typically do not contribute significantly are good candidates to neglect to reduce the number of high-quality ab initio calculations, as this allows us to reduce the stencil if some derivatives are not needed. Going to larger molecules, where the number of derivative terms increases due to the total number of sites, this is crucial to maintain computational feasibility.

Figure 3 shows the impact of the individual derivatives grouped by their order and combination of electronic, spatial, and charge dimensions as obtained from one million model estimations. This allows us to compare which interaction terms are relevant for the PECD parameters βi. We first notice that higher orders typically contribute substantially less: First order terms typically contribute 10%–35% to the prediction value, second order terms 5%–10%, and higher orders typically contribute less than 3%. This is reassuring, as this motivates the truncation of the series expansion after a few orders. For design applications, minor deviations from the true value are acceptable since the optimization can be performed using the model and then validated using a few accurate ab initio calculations.

FIG. 3.

Contribution of the different (mixed) partial derivative terms to the property values βi for up to fourth order. The variables with respect to which the property has been differentiated are grouped by electronic (E), spatial (R), and charge (Q) coordinates. Numbers indicate how many variables fall into the particular category, i.e., 2β1/∂Rx∂Qy is shown next to the label “R1 Q1.” The histograms are aggregated over all derivatives of identical labels. For readability, each histogram is normalized to the same peak value, and horizontal lines group terms by expansion order, as emphasized by common coloring. Data are shown for 106 random distortions chosen uniformly randomly over the search space. The bars show the 33rd to 67th percentile for each distribution. The tail of each distribution extends to the most extreme values and, therefore, represents the full domain of contributions for that group of derivatives.

FIG. 3.

Contribution of the different (mixed) partial derivative terms to the property values βi for up to fourth order. The variables with respect to which the property has been differentiated are grouped by electronic (E), spatial (R), and charge (Q) coordinates. Numbers indicate how many variables fall into the particular category, i.e., 2β1/∂Rx∂Qy is shown next to the label “R1 Q1.” The histograms are aggregated over all derivatives of identical labels. For readability, each histogram is normalized to the same peak value, and horizontal lines group terms by expansion order, as emphasized by common coloring. Data are shown for 106 random distortions chosen uniformly randomly over the search space. The bars show the 33rd to 67th percentile for each distribution. The tail of each distribution extends to the most extreme values and, therefore, represents the full domain of contributions for that group of derivatives.

Close modal

For both parameters, the electronic derivatives, i.e., the dependency on the photoelectron energy to first order, are by far the largest contribution. This means that future applications to molecules will surely require data from multiple photoelectron energy values. It is interesting to note that this term loses importance with higher orders quite quickly. Since there is only one input dimension related to the photoelectron energy, all even order contributions can either be positive or negative, as is evident from Fig. 3. Even though the distributions have a long tail, this tail is particularly thin. This likely means that the inclusion of third order energy derivatives is sufficient.

For the charge parameters, two orders cover almost all effects. It is interesting to see that for β2, all second order derivatives of the charges are positive. This means that the model geometry is a local minimum in β2. For β1, a highly asymmetric distribution in the second order contributions indicates the same effect for part of the charge dimensions. Higher order terms have negligible contributions. Therefore, second order derivatives are likely sufficient for molecular design.

For the spatial contributions, only first order terms are significant, which is particularly helpful for molecular design since this means that the response of βi due to the changed geometry of the design target is comparably easy to capture with a linear response. This is fortunate because there are three times as many spatial dimensions as there are charge dimensions to consider. Since obtaining the approximate minimum energy geometry for molecules is much cheaper than obtaining new PECD data via ab initio calculations, this offers a clear path toward a quick assessment of the molecular PECD parameters.

Since the mixed terms that couple different kinds of degrees of freedom may exhibit correlations, we analyzed their impact from ablation, i.e., dropping all groups of derivative contributions from the full set until the mean absolute error for a quantity reached 20% of the average quantity value. We expect this to be sufficiently accurate for iterative optimization and in line with the accuracy estimation from the first part of the results. We find that the second order spatial and charge degrees of freedom mix with the photoelectron energy but not among each other for both parameters. In the third order, only linear response, spatial, and charge parameters are coupled with the photoelectron energy. For β2, the linear response of the nuclear charges couples strongly with the second order photoelectron energy response.

Because each derivative coefficient requires one more ab initio calculation, it is desirable to define the minimal set of derivatives that is needed to obtain a quantum alchemy prediction with predictive power for β1. This is amplified in the case of treating real molecules with N sites to change, since there are 3N spatial and N charge degrees of freedom to consider. Each of the derivative terms has a different scaling with N: There is always a constant number of derivatives of order k with respect to the photoelectron energy, but, e.g., (3N)!/(3Nk)! many with respect to spatial degrees of freedom. Assuming a molecule to have ten atoms that could be altered, Fig. 4 shows the ablation results if individual terms are neglected in relation to the cost it would take to obtain these terms. There is a striking general trend that more expensive terms also contribute less—a strong argument for neglecting them to improve computational efficiency. Leaving out all terms in order of their computational cost until a mean error in the prediction of 20% is reached leaves 11 (mixed) derivatives, which for the system in this work require 29 ab initio calculations rather than the 85 done initially. For larger non-linear molecules with N atoms to change, there are 3N coefficients for all R1 terms: N for Q1, 3N for R1E1, N2 for Q2, N for Q1E1, 3N for R1E2, N for Q1E2, and 1 each for E0 to E4. Subtracting six coefficients to account for rotation and translation yields N2 + 12N − 1 calculations to determine all coefficients. Note that the results for the system in this work deviate from this formula since each of its atoms has only one spatial degree of freedom.

FIG. 4.

Error in the prediction value of β1 if a given derivative would be left out over the number of coefficients that would be needed for a molecule with ten atoms. Each derivative term is labeled and colored, as shown in Fig. 3. All terms that can be ignored until the error in the prediction reaches 20% are considered irrelevant. The inset shows a histogram of the actual prediction value and the approximated prediction value for 1.000 random points if only the terms marked relevant in the main figure are included. The inset shows that neglecting some of the most expensive derivative terms yields only little change in the overall prediction trend.

FIG. 4.

Error in the prediction value of β1 if a given derivative would be left out over the number of coefficients that would be needed for a molecule with ten atoms. Each derivative term is labeled and colored, as shown in Fig. 3. All terms that can be ignored until the error in the prediction reaches 20% are considered irrelevant. The inset shows a histogram of the actual prediction value and the approximated prediction value for 1.000 random points if only the terms marked relevant in the main figure are included. The inset shows that neglecting some of the most expensive derivative terms yields only little change in the overall prediction trend.

Close modal

This formula allows us to estimate the molecular size for which this perturbative approach is more efficient compared to brute force enumeration of the target search space. For energies, density matrix elements, orbital eigenvalues, and one-electron properties in real molecules, quantum alchemy typically20 converges for at least ΔZ = ±1. Here, this would allow three elements per atom, which can be varied during the search. For N changing sites, this allows for 3N targets to be estimated from the same reference perturbations. This search space clearly grows faster than the quadratic number of derivatives that need to be determined. With four atoms, the quantum alchemy approach becomes faster, and with nine atoms, it is already 100 times faster. It is worth noting that for larger systems, the overall electronic structure’s response to individual local changes becomes less pronounced. This would likely increase the convergence of the perturbation, which is in favor of the method.

In the present work, bond lengths and electron density distribution parameters of the model chiral system were optimized to achieve the desired values of the PECD. The screening of the variational space was performed by the Taylor series expansion algorithm using a relatively small set of ab initio computed data around the reference configuration. We were able to increase the dichroic parameter β1 by 3 to 4 times as compared to the reference geometry and, for some photoelectron energies, to achieve its value of around 26%. Keeping in mind that the dichroic parameter emerges on a typical scale of a few percent,1–3 the obtained huge value is rather unusual for randomly oriented targets, where sizable contributions from different molecular orientations (on the order of 20%–50%15–18) compensate each other almost completely, resulting in a much smaller effect. In addition, we were able to predict the parameter set at which the PECD changes its sign.

The results allowed us to dissect the contributions of different physical interactions of spatial and electronic degrees of freedom with the charge distributions, and that analysis offers clear information about which derivatives are likely to be of importance for molecular design to optimize the PECD parameters. We found that all expensive high-order contributions have only little impact on the overall prediction, and, therefore, the perturbative approach offers a computationally advantageous setting for systems with many atoms in the design space. This is particularly interesting as—unlike other electronic properties of molecules studied thus far—all degrees of freedom couple strongly for PECD parameters. This means that molecular design likely requires an iterative procedure where changes to the composition and geometry are suggested repeatedly in order to obtain large changes. The current work, however, showed that the one-step optimization already achieved a tremendous increase in the dichroic parameter.

In addition to the investigation of molecules, future work might explore whether theoretical limit expressions can be derived for βi, which could then be used in Padé approximants rather than Taylor series, or whether variable transformations allow for the disentanglement of the interplay of all degrees of freedom. Both routes would improve predictive accuracy, allowing further reduction of the computational cost. Automatic differentiation as applied to other ab initio calculations33–36 may be a promising alternative to finite difference evaluation, further reducing computational costs for design applications.

The supplementary material for this article includes the full set of the ab initio β1 and β2 parameters, computed at different initial sets of the variational parameters and used for the present optimization procedure.

A.N.A. and P.V.D. acknowledge the support from the Deutsche Forschungsgemeinschaft (DFG)—Project No. 328961117—SFB 1319 ELCH (Extreme light for sensing and driving molecular chirality, subproject C1).

The authors have no conflicts to disclose.

Guido F. von Rudorff: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Anton N. Artemyev: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Boris M. Lagutin: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Philipp V. Demekhin: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

2.
L.
Nahon
,
G. A.
Garcia
, and
I.
Powis
,
J. Electron Spectrosc. Relat. Phenom.
204
,
322
(
2015
).
3.
S.
Turchini
,
J. Phys.: Condens. Matter
29
,
503001
(
2017
).
5.
N.
Böwering
,
T.
Lischke
,
B.
Schmidtke
,
N.
Müller
,
T.
Khalil
, and
U.
Heinzmann
,
Phys. Rev. Lett.
86
,
1187
(
2001
).
6.
G. A.
Garcia
,
L.
Nahon
,
M.
Lebech
,
J. C.
Houver
,
D.
Dowek
, and
I.
Powis
,
J. Chem. Phys.
119
,
8781
(
2003
).
7.
C.
Lux
,
M.
Wollenhaupt
,
T.
Bolze
,
Q.
Liang
,
J.
Köhler
,
C.
Sarpe
, and
T.
Baumert
,
Angew. Chem., Int. Ed.
51
,
5001
(
2012
).
8.
C. S.
Lehmann
,
N. B.
Ram
,
I.
Powis
, and
M. H. M.
Janssen
,
J. Chem. Phys.
139
,
234307
(
2013
).
9.
S.
Beaulieu
,
A.
Ferré
,
R.
Géneaux
,
R.
Canonge
,
D.
Descamps
,
B.
Fabre
,
N.
Fedorov
,
F.
Légaré
,
S.
Petit
,
T.
Ruchon
,
V.
Blanchet
,
Y.
Mairesse
, and
B.
Pons
,
New J. Phys.
18
,
102002
(
2016
).
10.
11.
N.
Berova
,
K.
Nakanishi
, and
R. W.
Woody
,
Circular Dichroism: Principles and Applications
(
Wiley
,
New York
,
2009
).
12.
A.
Kastner
,
C.
Lux
,
T.
Ring
,
S.
Züllighoven
,
C.
Sarpe
,
A.
Senftleben
, and
T.
Baumert
,
ChemPhysChem
17
,
1119
(
2016
).
13.
R. E.
Goetz
,
C. P.
Koch
, and
L.
Greenman
,
Phys. Rev. Lett.
122
,
013204
(
2019
).
14.
R. E.
Goetz
,
C. P.
Koch
, and
L.
Greenman
,
J. Chem. Phys.
151
,
074106
(
2019
).
15.
M.
Tia
et al,
J. Phys. Chem. Lett.
8
,
2780
(
2017
).
16.
G.
Nalin
et al,
Phys. Chem. Chem. Phys.
23
,
17248
(
2021
).
17.
18.
19.
G. F.
von Rudorff
and
O. A.
von Lilienfeld
,
Phys. Rev. Res.
2
,
023220
(
2020
).
20.
G. F.
von Rudorff
,
J. Chem. Phys.
155
,
224103
(
2021
).
21.
E. A.
Eikey
,
A. M.
Maldonado
,
C. D.
Griego
,
G. F.
von Rudorff
, and
J. A.
Keith
,
J. Chem. Phys.
156
,
064106
(
2022
).
22.
A. N.
Artemyev
,
A. D.
Müller
,
D.
Hochstuhl
, and
P. V.
Demekhin
,
J. Chem. Phys.
142
,
244105
(
2015
).
23.
A. N.
Artemyev
,
E.
Kutscher
, and
P. V.
Demekhin
,
J. Chem. Phys.
156
,
031101
(
2022
).
24.
N. A.
Cherepkov
,
J. Phys. B: At. Mol. Phys.
14
,
2165
(
1981
).
N. A.
Cherepkov
,
Chem. Phys. Lett.
87
,
344
(
1982
).
25.
A.
Knie
,
M.
Ilchen
,
P.
Schmidt
,
P.
Reiß
,
C.
Ozga
,
B.
Kambs
,
A.
Hans
,
N.
Müglich
,
S. A.
Galitskiy
,
L.
Glaser
,
P.
Walter
,
J.
Viefhaus
,
A.
Ehresmann
, and
P. V.
Demekhin
,
Phys. Rev. A
90
,
013416
(
2014
).
26.
M.
Ilchen
,
G.
Hartmann
,
P.
Rupprecht
,
A. N.
Artemyev
,
R. N.
Coffee
,
Z.
Li
,
H.
Ohldag
,
H.
Ogasawara
,
T.
Osipov
,
D.
Ray
,
P.
Schmidt
,
T. J. A.
Wolf
,
A.
Ehresmann
,
S.
Moeller
,
A.
Knie
, and
P. V.
Demekhin
,
Phys. Rev. A
95
,
053423
(
2017
).
27.
A. D.
Müller
,
A. N.
Artemyev
, and
P. V.
Demekhin
,
J. Chem. Phys.
148
,
214307
(
2018
).
28.
A. D.
Müller
,
E.
Kutscher
,
A. N.
Artemyev
, and
P. V.
Demekhin
,
J. Chem. Phys.
152
,
044302
(
2020
).
29.
P. V.
Demekhin
,
A. N.
Artemyev
,
A.
Kastner
, and
T.
Baumert
,
Phys. Rev. Lett.
121
,
253201
(
2018
).
30.
E.
Kutscher
,
A. N.
Artemyev
, and
P. V.
Demekhin
,
Phys. Rev. A
107
,
013107
(
2023
).
31.
G. F.
von Rudorff
(
2024
). “
Nablachem/nablachem: 24.2
,” ,
2024
.
32.
M.
Baer
, maroba/findiff: Python package for numerical derivatives and partial differential equations in any number of dimensions,
2018
, https://github.com/maroba/findiff.
33.
M. F.
Kasim
,
S.
Lehtola
, and
S. M.
Vinko
,
J. Chem. Phys.
156
,
084801
(
2022
).
34.
A. S.
Abbott
,
B. Z.
Abbott
,
J. M.
Turney
, and
H. F.
Schaefer
III
,
J. Phys. Chem. Lett.
12
,
3232
(
2021
).
35.
X.
Zhang
and
G. K.-L.
Chan
,
J. Chem. Phys.
157
,
204801
(
2022
).
36.
T.
Tamayo-Mendoza
,
C.
Kreisbeck
,
R.
Lindh
, and
A.
Aspuru-Guzik
,
ACS Cent. Sci.
4
,
559
(
2018
).