Next Article in Journal
Synthesis and Self-Assembling Properties of Peracetylated β-1-Triazolyl Alkyl D-Glucosides and D-Galactosides
Next Article in Special Issue
Magnetic Aromaticity of Cycloporphyrin Nanorings
Previous Article in Journal
Dispersion of Micro Fibrillated Cellulose (MFC) in Poly(lactic acid) (PLA) from Lab-Scale to Semi-Industrial Processing Using Biobased Plasticizers as Dispersing Aids
Previous Article in Special Issue
Magnetic Shielding Study of Bonding and Aromaticity in Corannulene and Coronene
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Topological Analysis of Magnetically Induced Current Densities in Strong Magnetic Fields Using Stagnation Graphs

1
School of Chemistry, University of Nottingham, University Park, Nottingham NG7 2RD, UK
2
Hylleraas Centre for Quantum Molecular Sciences, Department of Chemistry, University of Oslo, P.O. Box 1033 Blindern, N-0315 Oslo, Norway
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Submission received: 6 August 2021 / Revised: 23 August 2021 / Accepted: 24 August 2021 / Published: 26 August 2021

Abstract

:
Stagnation graphs provide a useful tool to analyze the main topological features of the often complicated vector field associated with magnetically induced currents. Previously, these graphs have been constructed using response quantities appropriate for modest applied magnetic fields. We present an implementation capable of producing these graphs in arbitrarily strong magnetic fields, using current-density-functional theory. This enables us to study how the topology of the current vector field changes with the strength and orientation of the applied magnetic field. Applications to CH 4 , C 2 H 2 and C 2 H 4 are presented. In each case, we consider molecular geometries optimized in the presence of the magnetic field. The stagnation graphs reveal subtle changes to this vector field where the symmetry of the molecule remains constant. However, when the electronic state and symmetry of the corresponding equilibrium geometry changes with increasing field strength, the changes to the stagnation graph are extensive. We expect that the approach presented here will be helpful in interpreting changes in molecular structure and bonding in the strong-field regime.

1. Introduction

The topology and properties of magnetically induced currents have been widely studied using linear response techniques [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19] and a number of programs have been developed to both calculate the induced currents and analyze their main features [18,20,21,22,23]. An essential requirement is to ensure gauge–origin independence of the calculated currents; this has been achieved by a range of methods including: the individual gauge for localized orbitals (IGLO) method [24], the continuous set of gauge transformations (CSGT) approach [6,25,26], the continuous transformation of the gauge origin of the current density (CTOCD) method [8,10,14,15,27,28,29,30] and using London atomic orbitals (LAOs) [18,31] (also known as gauge-including atomic orbitals (GIAOs)). The physical current density is a rich source of chemical information [21,26,32,33,34,35,36,37]; its topology reflects the chemical structure of the molecule and its interaction with the applied magnetic field. Its features have been used as a criterion to assess the aromaticity of molecules [38], give insight into electron delocalization [39,40,41,42] and to probe hydrogen bond strengths [43,44,45]. The magnetically induced currents can also be directly related to magnetic response properties via the induced current susceptibility, a tensor describing the derivative of the induced current with respect to the applied magnetic field. Evaluating this response quantity at a zero magnetic field and using the Biot-Savart law leads directly to magnetic susceptibilities and NMR shielding constants. Furthermore, ring-current models have long been used to rationalize NMR chemical shifts in molecules [21,32,46,47,48,49].
The magnetically induced current is a relatively complicated vector field, and as such, tools for its analysis and interpretation of its main features in a simple manner are highly desirable. Approaches that analyze the induced currents by integration are well developed, see, for example, the GIMIC program [18,41]. Recently, we presented an implementation of these techniques in the context of the current-density-functional theory (CDFT), allowing applications to systems in arbitrarily strong magnetic fields [50]. Somewhat less attention has been given to topological approaches, which employ concepts from vector field analysis to analyze the topology of the vector field and provide insight into the magnetic behaviour of the system.
A few groups have pursued a topological analysis following early work on the mathematical characterization of stagnation points in vector fields by Reyn [51]. In particular, see the work of Keith and Bader [26], as well as works by Pelloni, Faglioni, Zanasi and Lazzeretti [52,53,54,55] for applications to molecular systems. These studies demonstrate how the location and classification of stagnation points (points in space at which the current density j has a magnitude of zero, i.e., | j | = 0 ) to produce stagnation graphs can distill the main features of the complicated current vector field into simpler plots that can be easily visualized without suffering from issues, such as occlusion, that often make the direct visualization of dense 3D vector fields challenging.
In the present work, we extend these techniques to analyze magnetically induced currents in strong magnetic fields at the CDFT level. This allows us to identify the main features and topology of induced currents as a function of the applied magnetic field strength and its direction relative to the molecular frame. In Section 2, we outline the necessary theoretical background to calculate the magnetically induced current densities in strong fields; detailing how the current density is determined in Section 2.1 and how its topological characteristics may be classified in Section 2.2. In Section 3, we describe the computational approach used to locate and classify stagnation points and construct stagnation graphs. The results are presented in Section 4 for applications to the CH 4 , C 2 H 2 and C 2 H 4 molecules at their ground state equilibrium geometries over a range of field strengths, obtained using a recently developed implementation of geometrical gradients for calculations using LAOs [56]. The changes in the stagnation plots with the applied field strength are carefully examined and analyzed in Section 5. Conclusions and directions for future work are discussed in Section 6.

2. Theoretical Background

2.1. Magnetically Induced Current Densities

In the presence of a uniform magnetic field B , the non-relativistic electronic Hamiltonian is given by
H ^ = H ^ 0 + 1 2 ( B × r O ) · p ^ + B · s ^ + 1 8 ( B × r O ) · ( B × r O ) ,
where H ^ 0 is the zero–field Hamiltonian, p ^ the momentum operator ( ı ), s ^ the spin operator and r O the position relative to some gauge–origin O . Since · B = 0 for a uniform magnetic field, the magnetic field may be represented by a vector potential A such that B = × A ; in the Coulomb gauge, this vector potential is defined to have a divergence of zero, · A = 0 . Therefore, for a uniform magnetic field, the vector potential depends on the gauge–origin as
A O ( r ) = 1 2 B × ( r O ) ,
and a change of the position of the gauge–origin O G is a gauge transformation,
A G ( r ) = A O ( r ) A O ( G ) · r .
This gauge–transformation corresponds to a unitary transform of the Hamiltonian,
H ^ = exp ( ı A O ( G ) · r ) H ^ exp ( ı A O ( G ) · r ) ,
the eigenfunctions of which, ψ , therefore undergo a compensating unitary transformation
ψ = exp ( ı A O ( G ) · r ) ψ ,
under which observables of the system remain gauge–origin invariant. This gauge–origin dependence of the wavefunction cannot be properly represented in a finite basis except by explicitly including the gauge–origin; this is the approach taken using LAOs [31], which comprise a standard Gaussian-type basis function φ a centred on R and multiplied by a field-dependent complex phase factor,
ω a ( r ) = φ a ( r ) exp ı 2 B × ( R O ) · r
which yields wavefunctions that are correct to first order in the magnetic field and properties that are gauge–origin invariant [57]. Utilizing an LAO basis, the effects of the magnetic field can be treated in a non-perturbative manner, allowing the behaviour of the systems in the magnetic fields of arbitrary strength to be examined [58,59].
The magnetically induced physical current density j is a continuous vector field in three dimensions and may be written as the sum of the diamagnetic current density j d and paramagnetic current density j p  [59,60],
j d ( r ) = A ( r ) σ ρ σ ( r )
j p ( r ) = ı 2 σ i ϕ i σ ( r ) ϕ i σ ( r ) ϕ i σ ( r ) ϕ i σ ( r )
where ρ σ is the σ -spin density, and ϕ i σ are the σ -spin occupied molecular orbitals.
Through the non-perturbative inclusion of magnetic field effects, the magnetically induced current density can be evaluated without the need for linear response calculations; in a basis of LAOs, the one–particle density matrix D a b σ computed at the Hartree–Fock (HF) or CDFT levels [59,60,61] may be used to evaluate the diamagnetic and paramagnetic current densities as
j d ( r ) = A ( r ) σ a b D a b σ ω a ( r ) ω b ( r )
j p ( r ) = ı 2 σ a b D a b σ ω a ( r ) ω b ( r ) ω a ( r ) ω b ( r ) ,
each of which are individually gauge–origin dependent; however, the physical current density j = j d + j p is gauge–origin invariant and can be computed at arbitrary field strengths.

2.2. Topological Characteristics

The magnetically induced current density j ( r ) is a three-dimensional vector field with a topological structure that may be characterized in terms of its singularities, otherwise referred to as stagnation points, at which | j | = 0 . The collection of these points for a system with magnetically induced currents is its stagnation graph, which describes the topological structure of the vector field [6,36,37,51,62]. The Cartesian components of the current density j α ( r ) at position r around a stagnation point r 0 may be described by the second–order Taylor expansion
j α ( r ) = β ( r β r 0 β ) j α r β r 0 + 1 2 β γ ( r β r 0 β ) ( r γ r 0 γ ) 2 j α r β r γ r 0 ,
in which the zeroth–order term is, by definition, zero. Taking only the first–order approximation, the current density in the region of the singular point can be described by the linear equation [54,63,64,65,66,67]:
j = J d , J α β = j α r β r 0 , d = r r 0 ,
where j is the current density at r , and J is the Jacobian matrix with elements J α β comprising the first derivative of j α with respect to r β .
It can be shown that the local behaviour of the current may be characterized by the eigenvalues η i of the Jacobian matrix [51]. The number of non-zero eigenvalues of J is denoted by the rank r, whilst the excess of eigenvalues with a positive real component over those with a negative real component is denoted the signature s – together the ordered pair ( r , s ) can be used to characterize the stagnation point [22,52,53,54,64,65]. Given that, at points where the current density is zero, the identity · j = 0 must be satisfied, the only possible ( r , s ) combinations for a three-dimensional vector field are ( 3 , ± 1 ) , ( 2 , 0 ) and ( 0 , 0 ) . In addition, a topological index i may be defined at a stagnation point where J has two non–zero eigenvalues as
i = + 1 η 1 = η 2 , η 3 = 0 1 η 1 = η 2 , η 3 = 0
The resulting classifications are summarized in Table 1 and will be used throughout this work [22,54,65,66,67].

3. Computational Methods

Previous works [22,65] have suggested using Newton–Raphson approaches to search for stagnation points, which occur at the nodes of the current density j . In particular, these optimization schemes only use first-order information to search the three-dimensional current-density vector field from a grid of arbitrarily chosen starting points [22]. Here, we present an alternative choice of the objective function in stagnation point searches in Section 3.1 before outlining an approach allowing for full second-order trust-region optimization, which benefits from quadratic convergence in the vicinity of stagnation points, in Section 3.2. Our approach to selecting an initial grid of starting points for the search and their subsequent refinement to produce clear stagnation graphs is detailed in Section 3.3.

3.1. Selecting an Appropriate Objective Function

The purpose of the stagnation point search is to locate the set of points { r } at which the objective function
| j ( r ) | = j x 2 ( r ) + j y 2 ( r ) + j z 2 ( r ) 1 2 .
is zero. Such points can be located by searching for stationary points in | j | and selecting those at which | j ( r ) | = 0 . Previous works have suggested optimizing | j ( r ) | directly via the Newton–Raphson approach, which requires only evaluation of the objective function and its first derivative [22,65]. However, it is clear that | j ( r ) | will exhibit cusps at the stagnation points. To demonstrate this, we consider an example of the ethene molecule oriented in the y z -plane with the two carbon nuclei equidistant from the origin along the z-axis in Figure 1. Plotting | j | along the line 2.5 y 2.5 bohr at z = 2.3 bohr, we expect to observe three stagnation points in line with the plots of Ref. [54] when a field of 0.1 B 0 is applied parallel to the C–C bond axis ( B 0 = ħ e 1 a 0 2 = 2.3505 × 10 5 T ). These points are clearly visible in Figure 1a; however, it can also be seen in Figure 1a that the expected cusps are present at the stagnation points along this line. The first derivative of Equation (14) may be readily evaluated as
| j | r α = j x | j | j x r α + j y | j | j y r α + j z | j | j z r α ,
and the singularities associated with the factors 1 / | j | at the stagnation points are clearly present in Figure 1b, with the total derivative of the objective function being discontinuous at these points. In practice, we find that, using this objective function, the cusps associated with the stagnation points can be approached to a sufficient proximity that the optimization can be terminated; however, the rate of convergence is somewhat hindered.
A full second-order approach, in which the issues associated with the singularities of | j ( r ) | are avoided, may be straightforwardly formulated by considering at alternative objective function, | j ( r ) | 2 . Figure 2 shows plots of | j | 2 and its first derivative in the same region plotted for | j | and its derivatives in Figure 1; it can be seen that this objective function is continuous, and both the gradient and the Hessian of this function are well defined at the stagnation points. This objective function, its gradient and its Hessian may then be evaluated as
| j | 2 = j x 2 + j y 2 + j z 2
| j | 2 r α = γ j γ r α j γ + j γ j γ r α
2 | j | 2 r α r β = γ 2 j γ r α r β j γ + 2 j γ r α j γ r β + j γ 2 j γ r α r β .
The gradient components are displayed in Figure 2b and are well behaved as expected.
Figure 2. (a) The value of | j | 2 ; (b) The Cartesian gradient components γ | j | 2 and its magnitude plotted from y = 2.5 to 2.5 bohr at z = 2.3 bohr in the ethene molecule in the y z -plane with a magnetic field of 0.1 B 0 along the z-axis.
Figure 2. (a) The value of | j | 2 ; (b) The Cartesian gradient components γ | j | 2 and its magnitude plotted from y = 2.5 to 2.5 bohr at z = 2.3 bohr in the ethene molecule in the y z -plane with a magnetic field of 0.1 B 0 along the z-axis.
Chemistry 03 00067 g002
Furthermore, the Hessian may be readily evaluated affording full second-order optimization at a modest cost. The required partial derivatives with respect to the position may be evaluated in terms of the LAOs and density matrices as
j α = ı 2 σ a b D a b σ ω a r α ω b ω a ω b r α A α σ a b D a b σ ω a ω b j α r β = ı 2 σ a b D a b σ 2 ω a r α r β ω b + ω a r α ω b r β ω a r β ω b r α ω a 2 ω b r α r β
A α σ a b D a b σ ω a r β ω b + ω a ω b r β 1 2 ε α γ β B γ σ a b D a b σ ω a ω b
2 j α r β r γ = ı 2 σ a b D a b σ 3 ω a r α r β r γ ω b + 2 ω a r α r β ω b r γ + 2 ω a r α r γ ω b r β + ω a r α 2 ω b r β r γ 2 ω a r β r γ ω b r α ω a r β 2 ω b r α r γ ω a r γ 2 ω b r α r β ω a 3 ω b r α r β r γ 1 2 ε α γ β B γ σ a b D a b σ ω a r γ ω b + ω a ω b r γ 1 2 ε α β γ B β σ a b D a b σ ω a r β ω b + ω a ω b r β A α σ a b D a b σ 2 ω a r β r γ ω b + ω a r β ω b r γ + ω a r γ ω b r β + ω a 2 ω b r β r γ
These derivatives have been implemented in the Quest program [23]. The correctness of each contribution was verified by finite-difference calculations with respect to r α .

3.2. Optimization Algorithm

To locate the stagnation points, | j | 2 was minimized with respect to r using a trust region approach. Using this method, a quadratic model of the objective function
m k ( r k + d ) = | j | k 2 + d T g k + 1 2 d T H k d
is constructed around each point visited in the optimization r k . Here, d defines the step taken from the point r k , and | j | k 2 , g k and H k are the objective function, its gradient and Hessian at r k , respectively,
| j | k 2 = | j ( r k ) | 2 g k = | j | 2 r α r k H k = 2 | j | 2 r α r β r k
To ensure progress in the optimization, the step d is determined by solving the trust-region subproblem
min d m k ( r k + d ) , subject to | | d | | 2 Δ k ,
where Δ k is the trust-radius. The step cannot exceed the trust region, in which the quadratic model is expected to be reliable. At each iteration, the accuracy of the step from the quadratic model is monitored using the ratio of the actual change in | j | 2 to that predicted by the quadratic model. If the step does not produce a sufficient decrease in | j | 2 , the trust-radius is reduced. Alternatively, if the model is particularly accurate (as would be the case if it was close to a stationary point for example), then the trust-radius may be increased; if the model is reasonably accurate, the trust-radius is kept the same. This approach is detailed in Algorithm 1, where the ratio used to control the trust-radius is denoted γ k . At each step, we solve the trust–region subproblem Equation (24) using the Steihaug-Toint truncated conjugate gradient algorithm [68].
Algorithm 1 Trust Region optimization
  • Δ 0 = 0.1
  • for k = 0 , 1 , 2 , do
  •     obtain d k by solving Equation (24)
  •      γ k = | j ( r k ) | 2 | j ( r k + d k ) | 2 m k ( r k ) m k ( r k + d k )
  •     if  γ k < 0.25  then
  •          Δ k + 1 = max ( 0.25 Δ k , 1.0 × 10 8 )
  •     else if  γ k > 0.75 and | | d k | | = Δ k  then
  •          Δ k + 1 = min ( 2 Δ k , 0.3 )
  •     else
  •          Δ k + 1 = Δ k
  •     end if
  •      r k + 1 = r k + d k
  • end for
In practice, we observe rapid convergence from a wide range of starting points, with quadratic convergence in the local region. The optimization is terminated when the following convergence criteria are satisfied: the maximum value of the gradient | j | 2 is ⩽2 × 10 6 au, its root-mean-square is ⩽ 10 6 au, and its Euclidian norm is ⩽ 10 8 au, the norm of the change in r is ⩽ 3 × 10 5 au and its maximum value ⩽ 6 × 10 4 au and the norm of the change in | j | 2 is ⩽ 10 8 au. This allows us to clearly distinguish the required stagnation points in the current-density vector field.

3.3. Initial Point Selection

A key aspect of making a stagnation point search computationally tractable is how the initial points for the optimizations are selected. Few details in this regard are discussed in previous studies of stagnation plots [22,65]. In the present work, the stagnation point search is carried out after a converged CDFT calculation from initial points defined using an atom-centred quadrature grid, with angular coordinates given by the eleventh degree Lebedev quadrature [69] and radial coordinates given by the LMG method [70] with threshold on the relative error of the radial integral of 10 2 . This grid is of the type used in DFT calculations but is much sparser than would be required for a full numerical integration of quantities, such as the electron density and its gradient. It does, however, retain a similar structure to the full DFT integration grid, with more points found close to the nuclei and less points as the distance from the nuclei increases. This structure provides a set of starting points that are well placed to resolve the details of often more complex stagnation lines in the vicinity of nuclei whilst also sampling those further away. In the present work, we initiate stagnation point searches from this initial grid of points.
Typically, plotting the converged points reveals the structure of the stagnation graph, but points may be relatively sparsely placed along the stagnation lines. In order to increase the number of stagnation points located, a second refinement stage is carried out. The path between stagnation points from the initial search that are separated by 1.5 bohr or less is divided into segments of approximately 0.1 bohr and a new initial point created at the midpoint of each segment. Once these points are generated, they are filtered to remove any points that are within 0.05 bohr of another point in the set. Optimizations are then carried out from these points to refine the stagnation graph. This strategy was effective at locating a larger number of singular points, particularly along stagnation lines, whilst minimizing the computational cost. Finally, since | j | 2 decays towards zero as we move away from the molecule, stagnation points in regions of negligible charge density ( ρ < 10 3 au) are discarded to leave only those of interest within the molecular volume.

4. Results

To test the efficiency of this new implementation and investigate how the current vector field topology and associated stagnation graph changes in the presence of strong magnetic fields, we study three small molecules: CH 4 , C 2 H 2 and C 2 H 4 . For each case, we consider fields of magnitude | B | = 0.0 0.2 B 0 and optimize the molecular geometries using the analytic gradient implementation of Ref. [56] at the cTPSS/u-aug-cc-pVTZ level. Here, the prefix u- indicates that the basis set is used in its uncontracted form; uncontracted basis sets are used to provide greater flexibility to describe the response of the wavefunction to the magnetic field. The spherical-harmonic form of these Gaussian basis functions are used throughout.
To aid the efficiency of the calculations density-fitting is used for all calculations, with the AutoAux auxiliary basis set. This autogenerated auxiliary basis set is generated following the approach outlined in Ref. [71] and provides a conservative auxiliary basis set constructed by considering the product space of the primary orbital basis. This choice helps to ensure that the results are consistent over a wide range of field strengths. Here, we note that the TZ quality basis set employed should be sufficiently accurate to describe systems in field strengths | B | 0.2 B 0 —see, for example, Refs. [50,72] for discussion of this point.
The stagnation plots are visualized along with the nuclear framework using the QMole tool, built using the Plotly python package [73], and are shown in the following subsections in static form; html files are provided as Supporting Information, which allow the reader to explore the plots in 3D using any modern web browser.

4.1. CH 4

The energy and main geometrical parameters for the equilibrium structure of CH 4 , optimized in magnetic fields of | B | = 0.05 B 0 , 0.10 B 0 , 0.15 B 0 and 0.20 B 0 aligned in the lowest energy orientation parallel to one of the C–H bonds, are presented in Table 2. It can be seen that, as expected for a closed-shell molecule, the energy exhibits diamagnetic behaviour and increases with field strength.
In the absence of a magnetic field, the equilibrium geometry of CH 4 has the familiar T d point group with all C–H bond lengths and H–C–H bond angles equal. Upon the application of a magnetic field, the point group symmetry of the molecule becomes that of the molecule and magnetic field combined; this usually lowers the symmetry of the system since, of the symmetry operations at zero field, only rotation axes parallel to the field, mirror planes perpendicular to the field and inversion centres remain [74]. In the case of CH 4 with a magnetic field applied parallel to one of the C–H bonds, the only symmetry element remaining is the three-fold rotation axis parallel to the field; hence, the point group is reduced to C 3 .
It can be seen in Table 2 that, with an applied magnetic field, the optimized geometry distorts away from the tetrahedral structure to a lower symmetry arrangement; the length of the C–H bond along the C 3 axis becomes distinct from that of the other bonds. Similarly, the H–C–H angles involving the axial H become distinct from those involving only non-axial H atoms. In Table 2, these quantities are reported in pairs, with the first value corresponding to the axial case and the second the non-axial case.
It can be seen that the axial C–H bond becomes compressed with increasing field strength, reducing in length by ∼ 0.5 pm at | B | = 0.20 B 0 relative to zero field. The trend for the non-axial bonds is less clear since their length remains approximately constant over this range of field strengths. In addition, the axial H–C–H angles become slightly more acute and the non-axial H–C–H angles become slightly more obtuse with increasing field strength, with the non-axial H–C–H angles being around 1.1 smaller than the axial H–C–H angles at | B | = 0.2 B 0 .
The CH 4 stagnation plots in the range | B | = 0.05 0.2 B 0 are shown in Figure 3. The stagnation plot at | B | = 0.05 B 0 closely resembles the plot presented in Ref. [54], with saddle lines coloured blue, para- and dia-tropic vortex lines coloured red and green, respectively, and branching points coloured purple. The latter are located close to the C 3 axis as expected. The similarity with Ref. [54] confirms the accuracy of the present implementation for locating the stagnation points and also that the classifications in Table 1 are correctly implemented.
In general, the structures of the stagnation plots display only minor changes with increasing field strength in this range. However, two trends with increasing field strength may be observed; the outer saddle and diatropic vortex lines dilate to have a slightly larger radius around the carbon atom, whilst simultaneously, the inner paratropic vortex and saddle lines contract towards the carbon atom with increasing field strength. The stagnation plots reflect the fact that, as the field strength increases, the induced current becomes more intense and compact around the carbon atom.

4.2. C 2 H 2

In Table 3, the energies and main geometrical parameters for the equilibrium structures of C 2 H 2 , optimized in magnetic fields of | B | = 0.05 B 0 , 0.10 B 0 , 0.15 B 0 and 0.20 B 0 aligned in the lowest energy orientation perpendicular to the C–C bond, are presented. Over the range of fields considered, the ground state of the molecule changes from that with M s = 0 to M s = 2 ; this occurs at a field strength of around 0.12 B 0 . Consistent with Ref. [75], the equilibrium geometry of C 2 H 2 in the lowest energy triplet state at zero field has a cis structure, with both hydrogen atoms on the same side of the C–C bond. However, at a field strength of around 0.10 B 0 , the trans structure in which the two hydrogen atoms are on either side of the C–C bond becomes lower in energy than the cis structure; hence, the M s = 2 state has a trans equilibrium geometry at | B | > 0.1 B 0 , where it is the ground state. The energies of the optimized M s = 0 , cis- M s = 2 and trans- M s = 2 states are summarized in Figure 4.
The equilibrium geometry of the M s = 0 state of C 2 H 2 remains linear whilst it is the ground state; however, the presence of the magnetic field leads to a reduction in symmetry from D h to C 2 h since only the two-fold rotation axis parallel to the field, the mirror plane perpendicular to the field and the inversion centre are retained. It can be seen in Table 3 that, as the field strength increases, both the C–C and C–H bonds contract slightly. In comparison, the trans- M s = 2 state has an electronic configuration in which an α -spin electron in a bonding orbital has been excited to an anti-bonding orbital and undergone a spin flip to a β -spin electron. In the stronger fields considered here, the favorable interaction of the unpaired β -electrons with the external magnetic field is sufficient to make the M s = 2 states lower in energy than the M s = 0 state, whilst in the stronger fields, the trans conformation is favored over the cis. Consistent with these changes in occupation, the C–C bond lengthens significantly, and the C–H bonds also lengthen but to a lesser extent. The H–C–C bond angle becomes 118.5 at 0.15 B 0 and becomes more acute with the increasing field. In the presence of the field, the trans- M s = 2 state adopts an orientation such that the field lies in the molecular plane, and the point group symmetry of the system is reduced from C 2 h in the absence of a field to C i .
The stagnation plots for C 2 H 2 are presented in Figure 5. The plot at | B | = 0.05 B 0 exhibits the same features as that obtained from response calculations in Ref. [76]. In addition to the classes of the stagnation point identified for CH 4 , isolated saddle nodes are visible in the plane containing the internuclear axis and perpendicular to the applied field. The stagnation graph at 0.10 B 0 has a very similar structure to that at 0.05 B 0 . The picture at 0.15 B 0 and 0.20 B 0 is entirely different since the stagnation graph of the ground state at this field strength contains no stagnation lines but only isolated paratropic vortices and saddle points. These features will be discussed further in Section 5.
To illustrate how the stagnation plots neatly summarize the topology of the magnetically induced current vector field, contour plots of | j | in the x z and y z planes of Figure 5b are presented in Figure 6. On the left, the darkest blue features in the x z -plane, representing the smallest | j | , show where the stagnation lines are located. The stagnation line perpendicular to the C–C bond and passing through its midpoint and those that form loops encircling the nuclei can both be seen in Figure 6a. On the right, the darkest blue features in the y z -plane capture the central diatropic vortex line displayed as a ring around the C–C bond midpoint in Figure 5. In addition, the points intersecting the bond axis are clearly represented, along with the four isolated saddle node stagnation points. This demonstrates that the stagnation plots accurately and succinctly capture the main features of the topology of the complicated vector field associated with the magnetically induced current.

4.3. C 2 H 4

The energy and main geometrical parameters for the equilibrium structure of C 2 H 4 , optimized in magnetic fields of | B | = 0.05 B 0 , 0.10 B 0 , 0.15 B 0 and 0.20 B 0 aligned in the lowest energy orientation parallel to the C–C bond, are presented in Table 4. In this orientation, the zero-field point group of D 2 h is reduced to C 2 h in a magnetic field since only the two-fold rotation axis along the C–C bond parallel to the field, the mirror plane perpendicular to the field and the centre of inversion remain. Consistent with this symmetry, all of the C–H bonds remain equivalent even with increasing field strength. Between | B | = 0.05 B 0 and 0.20 B 0 , the C–C and C–H bonds are compressed by ∼0.3 and ∼0.7 pm, respectively, whilst at the same time, the H–C–H bond angles become more acute, reducing by 5 over this range of field strengths.
The stagnation plots for C 2 H 4 are presented in Figure 7. As for CH 4 , the plot at | B | = 0.05 B 0 exhibits the same features as that obtained from response calculations in Ref. [54]. The structure of the stagnation graph remains similar for all the field strengths considered here; as such, only that at | B | = 0.05 B 0 is presented in Figure 7.

5. Discussion

Stagnation graphs of the kind presented in Figure 3, Figure 5 and Figure 7 provide convenient spatial descriptions of the magnetically induced current densities in molecules, with which various magnetic properties can be predicted. Analysis of the stagnation graphs obtained by response calculations for CH 4 , C 2 H 2 and C 2 H 4 has been presented in Refs. [54,76]. We now examine how this analysis changes over the range of fields considered in the present work.
In the diamagnetic CH 4 molecule, the current flow around the edges of the molecule in the low-density regions is diatropic and perpendicular to B . At the centre of these rings of diatropic current flow is the C 3 rotation axis of the molecule and along which lies a diatropic vortex line. Approaching regions of the molecule in which the charge density is greater, the structure of the magnetically induced current density becomes more complex with multiple individual circulations or toroidal vortices. At the centre of each lies a paratropic or diatropic vortex line for paratropic or diatropic vortices, respectively; these branch from the axial vortex line above the carbon atom and converge to the axial vortex line below the carbon atom such that they form closed loops. In the same region, saddle lines form a closed loop between the branching points on the axial vortex line and lie at the points of zero current flow between adjacent vortices.
As described in Section 4, it can be seen in Figure 3 that the structure of the stagnation graph in CH 4 remains generally similar with increasing magnetic field strength; however, the outer ring of saddle and diatropic vortex lines around the carbon atom dilate, whilst the inner ring of saddle and paratropic vortex lines contract with increasing field strength. This may be confirmed by considering the current densities at different field strengths; streamlines of the current in the x y plane perpendicular to the magnetic field and at a height of z = 0.2 bohr relative to the carbon atom at the origin in CH 4 at | B | = 0.1 B 0 and 0.2 B 0 are shown in Figure 8a,b, respectively. At the higher field strength, the magnitude of the current flow around the toroidal vortices becomes larger, resulting in the diatropic stagnation lines moving further from the nucleus. At the same time, the paratropic vortices are compressed towards the nucleus; hence the paratropic vortex lines moving towards the nucleus.
As presented in Section 4, the ground state of the C 2 H 2 molecule changes between | B | = 0.10 B 0 and 0.15 B 0 from M s = 0 to M s = 2 and the equilibrium geometry changes from a linear structure to adopt a trans conformation. The stagnation graph changes completely, with only a few isolated stagnation points remaining. The pattern of stagnation points at this geometry and field strength may be understood by considering the magnitude of the current density in the plane of the molecule; this is presented in Figure 9a, where the stagnation points can be clearly identified. Since the M s = 2 state of C 2 H 2 is not closed shell and has a different number of α and β electrons, the current density of each is not necessarily the same. This can be seen clearly in Figure 9b,c, depicting the norm of the α and β current densities, respectively. In the α case, a line of zero current density loops between the two carbon nuclei and around each, whilst in the β case, separate lines of zero current density appear to encircle each of the carbon nuclei and a line of zero current density at the midpoint of the C–C bond parallel to the magnetic field is visible.
The overall stagnation graph describes points where the magnitude of the total current density is zero; since this is a non-negative quantity by definition, at each stagnation point, the magnitudes of both the α and β current densities must vanish. The total stagnation graph, therefore, represents the intersection of the sets of stagnation points that would be obtained for the α and β spin currents independently. Since the topology of the α and β spin currents are significantly different, the intersection of their stagnation points results in a small number of points, which are visible in Figure 9a and located in Figure 5c,d.
For C 2 H 4 , the stagnation graph exhibits very little change between | B | = 0.05 B 0 and 0.20 B 0 . This may be rationalized by noting that the symmetry of the molecule does not change with increasing field strength in the range studied here. The relationship between a molecule’s symmetry elements, particularly mirror planes, and its stagnation graph has been discussed in detail in Refs. [52,53,54], where it is shown that the presence and position of stagnation lines can be determined from mirror planes. Whilst the equilibrium C–C and C–H bond lengths and the H–C–H bond angles change with increasing field strength, the symmetry of the molecule remains constant, and as such, the features present in the stagnation graph remain the same, notwithstanding distortions with increasing field strength similar to those in CH 4 .
Previous studies of stagnation graphs have mainly focused on small molecules with high symmetry and have used first-order methods to determine the location of stagnation points [22,52,53,54,65]. We expect that the second-order optimization approach for locating the stagnation points presented here should allow stagnation graphs to be mapped-out efficiently for more complex systems. As demonstrated with C 2 H 2 , this approach allows the changes in the stagnation graph to be examined as the symmetry of the molecule and its state changes, which is expected to become essential to study systems in stronger magnetic fields [56].

6. Conclusions

A second-order optimization method for mapping out the stagnation graphs of molecular systems has been presented. In this approach, stagnation points are located by minimizing | j | 2 , which, in contrast to | j | , has well-defined first and second derivatives, enabling the stagnation graph to be elucidated efficiently and robustly for general systems.
In contrast to previous work in this area [52,53,54,76], the magnetic field effects are here treated in a non-perturbative manner, allowing stagnation graphs to be computed at arbitrary field strengths and the effect of varying field strength on the characteristics of the stagnation graphs to be examined. Furthermore, the changes in the stagnation graph arising due to the effect of the magnetic field on the equilibrium geometry of the molecules has been accounted for by optimizing the molecular geometries at each field strength using the implementation of Ref. [56] before computing the stagnation graphs. This approach has been applied to study the stagnation graphs of three small molecules: CH 4 , C 2 H 2 and C 2 H 4 , which have previously been considered using response calculations [54,76], across a range of magnetic field strengths. In weak fields, the results obtained with the present approach are consistent with these earlier results, indicating the reliability of the implementation.
Upon increasing the field strength, we observe that the extent to which the stagnation graphs change depends strongly on the how the symmetry of the molecular structure and electronic state are affected by the magnetic field. In cases where the electronic state and molecular symmetry remain the same as in the absence of a field, only subtle changes to the stagnation graphs are observed, such as contractions and dilations of the radii of closed stagnation line loops. These changes can be explained by considering the magnitude of the current densities, which generally increase with increasing field strength, increasing the radii of toroidal vortices and increasing the distance between their corresponding vortical stagnation lines.
In cases where increasing the field strength results in a change in the electronic state and an accompanying change in the symmetry of the molecular geometry, much more extensive changes to the stagnation graph are observed. For example, the change in the ground state of C 2 H 2 from M s = 0 to M s = 2 at increasing field strength completely alters the molecular structure, symmetry and hence the stagnation graph. The structure of the stagnation graph was rationalized by considering the individual α and β spin current densities, with the stagnation points for the total current being the intersection of the stagnation points for the α and β spin currents individually.
We expect that the second-order approach for determining stagnation plots presented in this work will become a useful tool for understanding the electronic structure of more complex molecules in strong magnetic fields. The present implementation allows flexibility to study stagnation graphs for a wide range of uniform magnetic fields in a non-perturbative manner and to resolve their α and β -spin contributions. In the present work, only small molecules with up to two carbon atoms have been considered; in later work, we will apply this method to study current densities in more diverse molecular systems, such as homo- and heterocyclic aromatic molecules and aromatic systems, for which the stagnation graphs are known to show a wider range of features [66]. In future work, we will also consider the generalization of this approach to non-uniform magnetic fields, as described in Ref. [77]. As noted by Pelloni and Lazzeretti [76], the interaction between toroidal vortices and the gradient of non-uniform magnetic fields can be examined with the aid of stagnation graphs, and these may provide useful insight into these effects on nuclear magnetic shielding as well as more exotic properties, such as molecular anapole moments.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/chemistry3030067/s1, interactive html stagnation plots: CH4_B005.html, CH4_B010.html, CH4_B015.html, CH4_B020.html, C2H2_B005.html, C2H2_B010.html, C2H2_B015.html, C2H2_B020.html and C2H4_B005.html.

Author Contributions

Conceptualization, A.M.T. and T.J.P.I.; methodology, A.M.T., A.G. and T.J.P.I.; software, A.M.T., A.G. and T.J.P.I.; validation, A.G. and T.J.P.I.; formal analysis, A.M.T., A.G. and T.J.P.I.; investigation, A.M.T., A.G. and T.J.P.I.; resources, A.M.T.; data curation, A.G. and T.J.P.I.; writing—original draft preparation, A.M.T., A.G. and T.J.P.I.; writing—review and editing, A.M.T., A.G. and T.J.P.I.; visualization, T.J.P.I.; supervision, A.M.T.; project administration, A.M.T.; funding acquisition, A.M.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the European Research Council under H2020/ERC Consolidator Grant top DFT (Grant No. 772259).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All data supporting the conclusions in this manuscript are presented in the body of the document.

Acknowledgments

We would like to thank Paolo Lazzeretti for interesting discussion on this topic. We are grateful for access to the University of Nottingham’s Augusta HPC service.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Laws, E.A.; Stevens, R.M.; Lipscomb, W.N. Magnetic Properties of AlH and N2 from Coupled Hartree–Fock Theory. J. Chem. Phys. 1971, 54, 4269–4278. [Google Scholar] [CrossRef]
  2. Ditchfield, R. Self-consistent perturbation theory of diamagnetism. Mol. Phys. 1974, 27, 789–807. [Google Scholar] [CrossRef]
  3. Lazzeretti, P.; Rossi, E.; Zanasi, R. Magnetic properties and induced current density in acetylene. Int. J. Quantum Chem. 1984, 25, 1123–1134. [Google Scholar] [CrossRef]
  4. Wolinski, K.; Hinton, J.F.; Pulay, P. Efficient implementation of the gauge-independent atomic orbital method for NMR chemical shift calculations. J. Am. Chem. Soc. 1990, 112, 8251–8260. [Google Scholar] [CrossRef]
  5. Gauss, J. Calculation of NMR chemical shifts at second-order many-body perturbation theory using gauge-including atomic orbitals. Chem. Phys. Lett. 1992, 191, 614–620. [Google Scholar] [CrossRef]
  6. Keith, T.A.; Bader, R.F. Calculation of magnetic response properties using a continuous set of gauge transformations. Chem. Phys. Lett. 1993, 210, 223–231. [Google Scholar] [CrossRef]
  7. Ruud, K.; Helgaker, T.; Kobayashi, R.; Jørgensen, P.; Bak, K.L.; Jensen, H.J.A. Multiconfigurational self-consistent field calculations of nuclear shieldings using London atomic orbitals. J. Chem. Phys. 1994, 100, 8178–8185. [Google Scholar] [CrossRef] [Green Version]
  8. Coriani, S.; Lazzeretti, P.; Malagoli, M.; Zanasi, R. On CHF calculations of second-order magnetic properties using the method of continuous transformation of origin of the current density. Theor. Chim. Acta 1994, 89, 181–192. [Google Scholar] [CrossRef]
  9. Gauss, J.; Stanton, J.F. Coupled-cluster calculations of nuclear magnetic resonance chemical shifts. J. Chem. Phys. 1995, 103, 3561–3577. [Google Scholar] [CrossRef]
  10. Zanasi, R.; Lazzeretti, P.; Malagoli, M.; Piccinini, F. Molecular magnetic properties within continuous transformations of origin of the current density. J. Chem. Phys. 1995, 102, 7150–7157. [Google Scholar] [CrossRef]
  11. Lee, A.M.; Handy, N.C.; Colwell, S.M. The density functional calculation of nuclear shielding constants using London atomic orbitals. J. Chem. Phys. 1995, 103, 10095–10109. [Google Scholar] [CrossRef]
  12. Cheeseman, J.R.; Trucks, G.W.; Keith, T.A.; Frisch, M.J. A comparison of models for calculating nuclear magnetic resonance shielding tensors. J. Chem. Phys. 1996, 104, 5497–5509. [Google Scholar] [CrossRef]
  13. Rauhut, G.; Puyear, S.; Wolinski, K.; Pulay, P. Comparison of NMR Shieldings Calculated from Hartree-Fock and Density Functional Wave Functions Using Gauge-Including Atomic Orbitals. J. Phys. Chem. 1996, 100, 6310–6316. [Google Scholar] [CrossRef]
  14. Zanasi, R. Coupled Hartree–Fock calculations of molecular magnetic properties annihilating the transverse paramagnetic current density. J. Chem. Phys. 1996, 105, 1460–1469. [Google Scholar] [CrossRef]
  15. Ligabue, A.; Pincelli, U.; Lazzeretti, P.; Zanasi, R. Current Density Maps, Magnetizability, and Nuclear Magnetic Shielding Tensors for Anthracene, Phenanthrene, and Triphenylene. J. Am. Chem. Soc. 1999, 121, 5513–5518. [Google Scholar] [CrossRef]
  16. Jusélius, J.; Sundholm, D. Ab initio determination of the induced ring current in aromatic molecules. Phys. Chem. Chem. Phys. 1999, 1, 3429–3435. [Google Scholar] [CrossRef]
  17. Lazzeretti, P. Assessment of aromaticity via molecular response properties. Phys. Chem. Chem. Phys. 2004, 6, 217–223. [Google Scholar] [CrossRef]
  18. Jusélius, J.; Sundholm, D.; Gauss, J. Calculation of current densities using gauge-including atomic orbitals. J. Chem. Phys. 2004, 121, 3952–3963. [Google Scholar] [CrossRef] [PubMed]
  19. Lin, Y.C.; Jusélius, J.; Sundholm, D.; Gauss, J. Magnetically induced current densities in Al42- and Al44- species studied at the coupled-cluster level. J. Chem. Phys. 2005, 122, 214308. [Google Scholar] [CrossRef] [PubMed]
  20. Fliegl, H.; Taubert, S.; Lehtonen, O.; Sundholm, D. The gauge including magnetically induced current method. Phys. Chem. Chem. Phys. 2011, 13, 20500. [Google Scholar] [CrossRef] [PubMed]
  21. Sundholm, D.; Fliegl, H.; Berger, R.J. Calculations of magnetically induced current densities: Theory and applications. WIREs Comput. Mol. Sci. 2016, 6, 639–678. [Google Scholar] [CrossRef]
  22. Monaco, G.; Summa, F.F.; Zanasi, R. Program Package for the Calculation of Origin-Independent Electron Current Density and Derived Magnetic Properties in Molecular Systems. J. Chem. Inf. Model. 2020, 61, 270–283. [Google Scholar] [CrossRef]
  23. QUEST. A Rapid Development Platform for Quantum Electronic Structure Techniques. 2017. Available online: quest.codes (accessed on 23 August 2021).
  24. Kutzelnigg, W.; Wüllen, C.; Fleischer, U.; Franke, R.; Mourik, T. The IGLO method. Recent developments. In Nuclear Magnetic Shieldings and Molecular Structure; Springer: Dordrecht, The Netherlands, 1993; pp. 141–161. [Google Scholar] [CrossRef]
  25. Keith, T.; Bader, R. Calculation of magnetic response properties using atoms in molecules. Chem. Phys. Lett. 1992, 194, 1–8. [Google Scholar] [CrossRef]
  26. Keith, T.A.; Bader, R.F.W. Topological analysis of magnetically induced molecular current distributions. J. Chem. Phys. 1993, 99, 3669–3682. [Google Scholar] [CrossRef]
  27. Lazzeretti, P.; Malagoli, M.; Zanasi, R. Coupled Hartree–Fock calculations of origin-independent magnetic properties of benzene molecule. J. Chem. Phys. 1995, 102, 9619–9625. [Google Scholar] [CrossRef]
  28. Soncini, A.; Teale, A.M.; Helgaker, T.; Proft, F.D.; Tozer, D.J. Maps of current density using density-functional methods. J. Chem. Phys. 2008, 129, 074101. [Google Scholar] [CrossRef]
  29. Soncini, A. Charge and Spin Currents in Open-Shell Molecules: A Unified Description of NMR and EPR Observables. J. Chem. Theory Comput. 2007, 3, 2243–2257. [Google Scholar] [CrossRef] [PubMed]
  30. Lazzeretti, P. Methods of continuous translation of the origin of the current density revisited. In Marco Antonio Chaer Nascimento; Springer: Berlin/Heidelberg, Germany, 2012; pp. 103–115. [Google Scholar] [CrossRef]
  31. London, F. Théorie quantique des courants interatomiques dans les combinaisons aromatiques. J. Phys. Radium 1937, 8, 397–409. [Google Scholar] [CrossRef] [Green Version]
  32. Lazzeretti, P. Ring currents. Prog. Nucl. Magn. Reson. Spectrosc. 2000, 36, 1–88. [Google Scholar] [CrossRef]
  33. Pople, J. Molecular orbital theory of aromatic ring currents. Mol. Phys. 1958, 1, 175–180. [Google Scholar] [CrossRef]
  34. McWeeny, R. Ring currents and proton magnetic resonance in aromatic molecules. Mol. Phys. 1958, 1, 311–321. [Google Scholar] [CrossRef]
  35. Hegstrom, R.A.; Lipscomb, W.N. Magnetic Properties of the BH Molecule. J. Chem. Phys. 1966, 45, 2378–2383. [Google Scholar] [CrossRef]
  36. Gomes, J.A.N.F. Topological elements of the magnetically induced orbital current densities. J. Chem. Phys. 1983, 78, 4585–4591. [Google Scholar] [CrossRef]
  37. Gomes, J.A.N.F. Topology of the electronic current density in molecules. Phys. Rev. A 1983, 28, 559–566. [Google Scholar] [CrossRef]
  38. Fliegl, H.; Sundholm, D.; Taubert, S.; Jusélius, J.; Klopper, W. Magnetically Induced Current Densities in Aromatic, Antiaromatic, Homoaromatic, and Nonaromatic Hydrocarbons. J. Phys. Chem. A 2009, 113, 8668–8676. [Google Scholar] [CrossRef]
  39. Taubert, S.; Sundholm, D.; Pichierri, F. Magnetically Induced Currents in Bianthraquinodimethane-Stabilized Möbius and Hückel [16]Annulenes. J. Org. Chem. 2009, 74, 6495–6502. [Google Scholar] [CrossRef]
  40. Fliegl, H.; Sundholm, D.; Taubert, S.; Pichierri, F. Aromatic Pathways in Twisted Hexaphyrins. J. Phys. Chem. A 2010, 114, 7153–7161. [Google Scholar] [CrossRef]
  41. Fliegl, H.; Sundholm, D.; Pichierri, F. Aromatic pathways in mono- and bisphosphorous singly Möbius twisted [28] and [30]hexaphyrins. Phys. Chem. Chem. Phys. 2011, 13, 20659. [Google Scholar] [CrossRef]
  42. Dickens, T.K.; Mallion, R.B.; Radenković, S. Assessing the Extent of π-Electron Delocalization in Naphtho-Annelated Fluoranthenes by Means of Topological Ring-Currents. J. Phys. Chem. A 2019, 123, 1445–1450. [Google Scholar] [CrossRef] [PubMed]
  43. Monaco, G.; Porta, P.D.; Jabłoński, M.; Zanasi, R. Topology of the magnetically induced current density and proton magnetic shielding in hydrogen bonded systems. Phys. Chem. Chem. Phys. 2015, 17, 5966–5972. [Google Scholar] [CrossRef] [PubMed]
  44. Fliegl, H.; Lehtonen, O.; Sundholm, D.; Kaila, V.R.I. Hydrogen-bond strengths by magnetically induced currents. Phys. Chem. Chem. Phys. 2011, 13, 434–437. [Google Scholar] [CrossRef]
  45. Porta, P.D.; Zanasi, R.; Monaco, G. Hydrogen-hydrogen bonding: The current density perspective. J. Comput. Chem. 2015, 36, 707–716. [Google Scholar] [CrossRef]
  46. Cao, J.; London, G.; Dumele, O.; von Wantoch Rekowski, M.; Trapp, N.; Ruhlmann, L.; Boudon, C.; Stanger, A.; Diederich, F. The Impact of Antiaromatic Subunits in [4n+2] π-Systems: Bispentalenes with [4n+2] π-Electron Perimeters and Antiaromatic Character. J. Am. Chem. Soc. 2015, 137, 7178–7188. [Google Scholar] [CrossRef] [PubMed]
  47. Baryshnikov, G.V.; Karaush, N.N.; Valiev, R.R.; Minaev, B.F. Aromaticity of the completely annelated tetraphenylenes: NICS and GIMIC characterization. J. Mol. Model. 2015, 21. [Google Scholar] [CrossRef] [PubMed]
  48. Baryshnikov, G.V.; Valiev, R.R.; Karaush, N.N.; Sundholm, D.; Minaev, B.F. Aromaticity of the doubly charged [8]circulenes. Phys. Chem. Chem. Phys. 2016, 18, 8980–8992. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Haigh, C.; Mallion, R. Ring current theories in nuclear magnetic resonance. Prog. Nucl. Magn. Reson. Spectrosc. 1979, 13, 303–344. [Google Scholar] [CrossRef]
  50. Irons, T.J.P.; Spence, L.; David, G.; Speake, B.T.; Helgaker, T.; Teale, A.M. Analyzing Magnetically Induced Currents in Molecular Systems Using Current-Density-Functional Theory. J. Phys. Chem. A 2020, 124, 1321–1333. [Google Scholar] [CrossRef]
  51. Reyn, J.W. Classification and description of the singular points of a system of three linear differential equations. Z. Angew. Math. Phys. 1964, 15, 540–557. [Google Scholar] [CrossRef]
  52. Pelloni, S.; Faglioni, F.; Zanasi, R.; Lazzeretti, P. Topology of magnetic-field-induced current-density field in diatropic monocyclic molecules. Phys. Rev. A 2006, 74, 012506. [Google Scholar] [CrossRef] [Green Version]
  53. Pelloni, S.; Lazzeretti, P.; Zanasi, R. Topological models of magnetic field induced current density field in small molecules. Theor. Chem. Acc. 2009, 123, 353–364. [Google Scholar] [CrossRef]
  54. Pelloni, S.; Lazzeretti, P. Stagnation graphs and topological models of magnetic-field induced electron current density for some small molecules in connection with their magnetic symmetry. Int. J. Quantum Chem. 2010, 111, 356–367. [Google Scholar] [CrossRef]
  55. Lazzeretti, P. Topological definition of ring currents. Phys. Chem. Chem. Phys. 2016, 18, 11765–11771. [Google Scholar] [CrossRef] [PubMed]
  56. Irons, T.J.P.; David, G.; Teale, A.M. Optimizing Molecular Geometries in Strong Magnetic Fields. J. Chem. Theory Comput. 2021, 17, 2166–2185. [Google Scholar] [CrossRef] [PubMed]
  57. Ditchfield, R. On molecular orbital theories of NMR chemical shifts. Chem. Phys. Lett. 1972, 15, 203–206. [Google Scholar] [CrossRef]
  58. Tellgren, E.I.; Soncini, A.; Helgaker, T. Nonperturbative ab initio calculations in strong magnetic fields using London orbitals. J. Chem. Phys. 2008, 129, 154114. [Google Scholar] [CrossRef]
  59. Tellgren, E.I.; Teale, A.M.; Furness, J.W.; Lange, K.K.; Ekström, U.; Helgaker, T. Non-perturbative calculation of molecular magnetic properties within current-density functional theory. J. Chem. Phys. 2014, 140, 034101. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  60. Vignale, G.; Rasolt, M. Density-functional theory in strong magnetic fields. Phys. Rev. Lett. 1987, 59, 2360–2363. [Google Scholar] [CrossRef]
  61. Furness, J.W.; Verbeke, J.; Tellgren, E.I.; Stopkowicz, S.; Ekström, U.; Helgaker, T.; Teale, A.M. Current Density Functional Theory Using Meta-Generalized Gradient Exchange-Correlation Functionals. J. Chem. Theory Comput. 2015, 11, 4169–4181. [Google Scholar] [CrossRef] [PubMed]
  62. Lazzeretti, P.; Rossi, E.; Zanasi, R. Singularities of magnetic-field induced electron current density: A study of the ethylene molecule. Int. J. Quantum Chem. 1984, 25, 929–940. [Google Scholar] [CrossRef]
  63. Lazzeretti, P. Continuity equations for electron charge densities and current densities induced in molecules by electric and magnetic fields. J. Chem. Phys. 2019, 151, 114108. [Google Scholar] [CrossRef]
  64. Lazzeretti, P. Stagnation graphs and separatrices of local bifurcations in velocity and current density planar vector fields. Rend. Lincei. Sci. Fis. Nat. 2019, 30, 515–535. [Google Scholar] [CrossRef]
  65. Monaco, G.; Zanasi, R. Magnetically Induced Current Density Spatial Domains. J. Phys. Chem. A 2019, 123, 1558–1569. [Google Scholar] [CrossRef]
  66. Pelloni, S.; Lazzeretti, P. Magnetotropicity of five-membered heterocyclic molecules. Theor. Chem. Acc. 2006, 117, 903–913. [Google Scholar] [CrossRef]
  67. Pelloni, S.; Lazzeretti, P.; Zanasi, R. Induced Orbital Paramagnetism and Paratropism in Closed-Shell Molecules. J. Phys. Chem. A 2009, 113, 14465–14479. [Google Scholar] [CrossRef] [PubMed]
  68. Steihaug, T. The Conjugate Gradient Method and Trust Regions in Large Scale Optimization. SIAM J. Sci. Comput. 1983, 20, 626–637. [Google Scholar] [CrossRef] [Green Version]
  69. Lebedev, V. Quadratures on a sphere. USSR Comput. Math. Math. Phys. 1976, 16, 10–24. [Google Scholar] [CrossRef]
  70. Lindh, R.; Malmqvist, P.Å.; Gagliardi, L. Molecular integrals by numerical quadrature. I. Radial integration. Theor. Chem. Acc. 2001, 106, 178–187. [Google Scholar] [CrossRef] [Green Version]
  71. Stoychev, G.L.; Auer, A.A.; Neese, F. Automatic Generation of Auxiliary Basis Sets. J. Chem. Theory Comput. 2017, 13, 554–562. [Google Scholar] [CrossRef]
  72. Lehtola, S.; Dimitrova, M.; Sundholm, D. Fully numerical electronic structure calculations on diatomic molecules in weak to strong magnetic fields. Mol. Phys. 2019, 118, e1597989. [Google Scholar] [CrossRef] [Green Version]
  73. Plotly. Plotly Python Graphing Library. 2021. Available online: plotly.com/python (accessed on 23 August 2021).
  74. Ceulemans, A.J. Group Theory Applied to Chemistry; Springer: Dordrecht, The Netherlands, 2013. [Google Scholar] [CrossRef]
  75. Wetmore, R.W.; Schaefer, H.F. Triplet electronic states of acetylene: Cis and trans structures and energetics. J. Chem. Phys. 1978, 69, 1648–1654. [Google Scholar] [CrossRef]
  76. Pelloni, S.; Lazzeretti, P. Ring current models for acetylene and ethylene molecules. Chem. Phys. 2009, 356, 153–163. [Google Scholar] [CrossRef]
  77. Tellgren, E.I.; Fliegl, H. Non-perturbative treatment of molecules in linear magnetic fields: Calculation of anapole susceptibilities. J. Chem. Phys. 2013, 139, 164118. [Google Scholar] [CrossRef] [PubMed]
Figure 1. (a) The value of | j | ; (b) The Cartesian gradient components γ | j | and its magnitude plotted from y = 2.5 to 2.5 bohr at z = 2.3 bohr in the ethene molecule in the y z -plane with a magnetic field of 0.1 B 0 along the z-axis.
Figure 1. (a) The value of | j | ; (b) The Cartesian gradient components γ | j | and its magnitude plotted from y = 2.5 to 2.5 bohr at z = 2.3 bohr in the ethene molecule in the y z -plane with a magnetic field of 0.1 B 0 along the z-axis.
Chemistry 03 00067 g001
Figure 3. The stagnation graphs of the CH 4 molecule in its equilibrium geometry with a magnetic field of (a) 0.05 B 0 , (b) 0.10 B 0 , (c) 0.15 B 0 and (d) 0.20 B 0 along the z-axis, arranged from left to right, respectively. The interactive version of these figures may be found in CH4_B005.html, CH4_B010.html, CH4_B015.html and CH4_B020.html of the Supporting Information, respectively.
Figure 3. The stagnation graphs of the CH 4 molecule in its equilibrium geometry with a magnetic field of (a) 0.05 B 0 , (b) 0.10 B 0 , (c) 0.15 B 0 and (d) 0.20 B 0 along the z-axis, arranged from left to right, respectively. The interactive version of these figures may be found in CH4_B005.html, CH4_B010.html, CH4_B015.html and CH4_B020.html of the Supporting Information, respectively.
Chemistry 03 00067 g003
Figure 4. The energies of the optimized M s = 0 , cis- M s = 2 and trans- M s = 2 states of the C 2 H 2 molecule with respect to the strength of magnetic field applied along the x-axis.
Figure 4. The energies of the optimized M s = 0 , cis- M s = 2 and trans- M s = 2 states of the C 2 H 2 molecule with respect to the strength of magnetic field applied along the x-axis.
Chemistry 03 00067 g004
Figure 5. The stagnation graphs of the C 2 H 2 molecule in its equilibrium geometry with a magnetic field of (a) 0.05 B 0 , (b) 0.10 B 0 , (c) 0.15 B 0 and (d) 0.20 B 0 along the x-axis. The interactive version of these figures may be found in C2H2_B005.html, C2H2_B010.html, C2H2_B015.html and C2H2_B020.html of the Supporting Information, respectively.
Figure 5. The stagnation graphs of the C 2 H 2 molecule in its equilibrium geometry with a magnetic field of (a) 0.05 B 0 , (b) 0.10 B 0 , (c) 0.15 B 0 and (d) 0.20 B 0 along the x-axis. The interactive version of these figures may be found in C2H2_B005.html, C2H2_B010.html, C2H2_B015.html and C2H2_B020.html of the Supporting Information, respectively.
Chemistry 03 00067 g005
Figure 6. Contour plots of the norm of the current density | j | in the (a) x z plane (left) and (b) y z plane (right) of the C 2 H 2 molecule at its equilibrium geometry with a magnetic field of 0.10 B 0 along the x-axis.
Figure 6. Contour plots of the norm of the current density | j | in the (a) x z plane (left) and (b) y z plane (right) of the C 2 H 2 molecule at its equilibrium geometry with a magnetic field of 0.10 B 0 along the x-axis.
Chemistry 03 00067 g006
Figure 7. The stagnation graph of the C 2 H 4 molecule in its equilibrium geometry with a magnetic field of 0.05 B 0 along the z-axis. The interactive version of this figure may be found in C2H4_B005.html of the Supporting Information.
Figure 7. The stagnation graph of the C 2 H 4 molecule in its equilibrium geometry with a magnetic field of 0.05 B 0 along the z-axis. The interactive version of this figure may be found in C2H4_B005.html of the Supporting Information.
Chemistry 03 00067 g007
Figure 8. Streamline plots of the current density j , with contours of its norm, in the x y plane of the CH 4 molecule in a magnetic field of (a) 0.1 B 0 and (b) 0.2 B 0 along the z-axis. Paratropic vortex lines are shown as red circles, diatropic vortex lines as green circles and saddle lines as blue circles.
Figure 8. Streamline plots of the current density j , with contours of its norm, in the x y plane of the CH 4 molecule in a magnetic field of (a) 0.1 B 0 and (b) 0.2 B 0 along the z-axis. Paratropic vortex lines are shown as red circles, diatropic vortex lines as green circles and saddle lines as blue circles.
Chemistry 03 00067 g008
Figure 9. Contour plots of the norm of the (a) current density | j | , (b) spin α current density | j α | and (c) spin β current density | j β | in the x z plane of the C 2 H 2 molecule at its equilibrium geometry with a magnetic field of 0.20 B 0 along the x-axis. In (a), paratropic vortex lines are shown as red circles and saddle lines as blue circles.
Figure 9. Contour plots of the norm of the (a) current density | j | , (b) spin α current density | j α | and (c) spin β current density | j β | in the x z plane of the C 2 H 2 molecule at its equilibrium geometry with a magnetic field of 0.20 B 0 along the x-axis. In (a), paratropic vortex lines are shown as red circles and saddle lines as blue circles.
Chemistry 03 00067 g009
Table 1. The classification of stagnation points by their topological characteristics.
Table 1. The classification of stagnation points by their topological characteristics.
( r , s ) i η i ConditionsMagnetotropicityClassification
( 3 , ± 1 ) η 1 , η 2 , η 3 R isolated singularity: saddle
( 3 , ± 1 ) η 3 R , η 1 = η 2 B · ( × j ) > 0 isolated singularity: paramagnetic focus
( 3 , ± 1 ) η 3 R , η 1 = η 2 B · ( × j ) < 0 isolated singularity: diamagnetic focus
( 2 , 0 ) 1 η 1 , η 2 R , η 3 = 0 stagnation line: saddle
( 2 , 0 ) + 1 Re ( η 1 ) , Re ( η 2 ) = 0 , η 3 = 0 B · ( × j ) > 0 stagnation line: paramagnetic vortex
( 2 , 0 ) + 1 Re ( η 1 ) , Re ( η 2 ) = 0 , η 3 = 0 B · ( × j ) < 0 stagnation line: diamagnetic vortex
( 0 , 0 ) η 1 , η 2 , η 3 = 0 branching point
Table 2. Equilibrium geometries of CH 4 at the magnetic field strengths considered in this work, optimized with the cTPSS functional and in the u-aug-cc-pVTZ basis. In all cases, the magnetic field is aligned parallel to one of the C–H bonds.
Table 2. Equilibrium geometries of CH 4 at the magnetic field strengths considered in this work, optimized with the cTPSS functional and in the u-aug-cc-pVTZ basis. In all cases, the magnetic field is aligned parallel to one of the C–H bonds.
| B | / B 0 Energy / E h R C H / bohr∢ H–C–H / DegreePoint Group
0.00−40.5415542.0638109.5T d
0.05−40.5366082.0631, 2.0632109.4, 109.5C 3
0.10−40.5221242.0614, 2.0615109.3, 109.6C 3
0.15−40.4991292.0585, 2.0604109.2, 109.8C 3
0.20−40.4692432.0543, 2.0620108.9, 110.0C 3
Table 3. Equilibrium geometries of C 2 H 2 at the magnetic field strengths considered in this work, optimized with the cTPSS functional and in the u-aug-cc-pVTZ basis. In all cases, the magnetic field is aligned perpendicular to the C–C bond.
Table 3. Equilibrium geometries of C 2 H 2 at the magnetic field strengths considered in this work, optimized with the cTPSS functional and in the u-aug-cc-pVTZ basis. In all cases, the magnetic field is aligned perpendicular to the C–C bond.
| B | / B 0 Energy / E h R C C / bohrR C H / bohr∢ H–C–C / DegreePoint Group
0.00−77.3744202.27272.0136180.0D h
0.05−77.3688532.27222.0123180.0C 2 h
0.10−77.3522242.27032.0085180.0C 2 h
0.15−77.3782852.59742.0679118.5C i
0.20−77.4217422.61062.0658115.8C i
Table 4. Equilibrium geometries of C 2 H 4 at the magnetic field strengths considered in this work, optimized with the cTPSS functional and in the u-aug-cc-pVTZ basis. In all cases, the magnetic field is aligned parallel to the C–C bond.
Table 4. Equilibrium geometries of C 2 H 4 at the magnetic field strengths considered in this work, optimized with the cTPSS functional and in the u-aug-cc-pVTZ basis. In all cases, the magnetic field is aligned parallel to the C–C bond.
| B | / B 0 Energy / E h R C C / bohrR C H / bohr∢ H–C–H / DegreePoint Group
0.00−78.6338152.51502.0532116.5D 2 h
0.05−78.6289542.51462.0523116.2C 2 h
0.10−78.6142442.51332.0497115.2C 2 h
0.15−78.5894252.51102.0457113.7C 2 h
0.20−78.5543912.50842.0407111.5C 2 h
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Irons, T.J.P.; Garner, A.; Teale, A.M. Topological Analysis of Magnetically Induced Current Densities in Strong Magnetic Fields Using Stagnation Graphs. Chemistry 2021, 3, 916-934. https://0-doi-org.brum.beds.ac.uk/10.3390/chemistry3030067

AMA Style

Irons TJP, Garner A, Teale AM. Topological Analysis of Magnetically Induced Current Densities in Strong Magnetic Fields Using Stagnation Graphs. Chemistry. 2021; 3(3):916-934. https://0-doi-org.brum.beds.ac.uk/10.3390/chemistry3030067

Chicago/Turabian Style

Irons, Tom J. P., Adam Garner, and Andrew M. Teale. 2021. "Topological Analysis of Magnetically Induced Current Densities in Strong Magnetic Fields Using Stagnation Graphs" Chemistry 3, no. 3: 916-934. https://0-doi-org.brum.beds.ac.uk/10.3390/chemistry3030067

Article Metrics

Back to TopTop