New LCAO-MO Method for Molecules in The Arbitrary Magnetic Field

Abstract

The Hartree-Fock-Slater (HFS) equation including magnetic interaction was solved full-numerically by the LCAO-MO method. Firstly, the AOs of each atom in a molecule were obtained by solving the HFS equation dedicated to the atom. This equation includes the same magnetic interaction as the equation of the molecule and virtual "well potential" to adjust the spreads of AOs. In this method, the Fock operator for the atom is expressed as a numerical matrix and the equation can be solved by a mathematical solver of the numerical eigenvalue problem. Then an MO is obtained by using these AOs as bases of the LCAO-MO method. AOs naturally have suitable complex phases for the MO to satisfy gauge invariance. To evaluate the new method, magnetic shielding constants σ were calculated for various small molecules comprising 2-3 atoms of second-row element and hydrogen atoms. The "well potential" of each atom was adjusted systematically and calculated σ showed fairly good agreement with experimental values.

1 Introduction

The LCAO-MO method is generally used to obtain MOs of various molecules. But MOs expressed in complex numbers, like MOs in magnetic fields, cause difficulty, because the functional space of such MO is too vast to express as a linear combination of a limited number of ordinary basis functions. In the case of the uniform magnetic field, this difficulty can be avoided by using London AO [1], which has a predetermined complex phase and only the real part of the AO needs to be calculated [2]. However, this method can't be applied in the case of an arbitrary magnetic field.

The author has proposed a new method where the HFS equation [3], which is prepared for each atom and includes magnetic interaction, is solved numerically. Obtained AOs are used as bases, which include complex phases corresponding to magnetic fields, so their linear combination is expected to express MO properly not only for uniform but also for arbitrary magnetic fields. But calculated AOs tend to spread too widely as bases. For such a case, the DV-Xα method [4], which solves the HFS equation numerically, prepared "well potential" [5] to adjust the spreads of AOs. Unfortunately, it contains adjustable parameters, the author has proposed a new method to optimize such parameters.

2 Method

The base equation of this method for a molecule is

  

$\begin{array}{c}(-\frac{1}{2}{\nabla }^2-\sum _A\frac{Z_A}{r_A}+\int \frac{{\rho }_2}{r_{12}}d{v}_2-X\alpha \\ +\frac{i}{C}A\cdot \nabla +\frac{1}{2}{\left(\frac{1}{C}\right)}^2{A}^2)\psi =\epsilon \psi \end{array}$(1)

where Xα designates Slater's Xα potential [3] to approximate the exchange interaction of Hartree Fock equation and A is vector potential. $\left(B=\nabla \times A\right)$ The Equation (1) is expressed in the atomic units [6], where m_e, e, $\hslash$, and k_e$\left(=1/4\pi {ϵ}_0\right)$ are defined as fundamental units, and their numerical values are set to be 1.0. C=137.036 in atomic units. The last 2 terms in parenthesis of (1) express magnetic interaction. MO ψ is expressed as a linear combination of AO χr which is a solution of the similar equation as (1) but its Z_A and ρ₂ include only a specific atom's contribution.

To treat arbitrary A, the equation is solved fully numerically. Sample points r_i (i = 1,2,…,n) are distributed densely near the nucleus and sparsely in the peripheral area for each atom. A set of χ_ri on each r_i constructs a state vector of the atom.

  

$\begin{array}{c}\left[\begin{array}{cccc}{F}_{11}& F_{12}& \cdots & F_{1n}\\ F_{21}& F_{22}& \cdots & F_{2n}\\ \cdots & \cdots & \ddots & ⋮\\ F_{n1}& F_{n2}& \cdots & F_{nn}\end{array}\right]\left[\begin{array}{c}{\chi }_{r1}\\ {\chi }_{r2}\\ ⋮\\ {\chi }_{rn}\end{array}\right]={\epsilon }_r\left[\begin{array}{c}{\chi }_{r1}\\ {\chi }_{r2}\\ ⋮\\ {\chi }_{rn}\end{array}\right]\end{array}$(2)

The base equation for the atom has the form (2) where [χ_ri] is a state vector and [F_ij] is a Fock operator. [F_ij] is the sum of [P] corresponding to each physical quantity P in the parenthesis of (1). Each matrix is constructed as follows.

A matrix that isn't related to ∇ is a diagonal matrix, whose diagonal element is P_i the value of P on r_i.

A matrix form of $[\nabla ] \left(=\left[\partial /\partial x\right],\left[\partial /\partial y\right],\left[\partial /\partial z\right]\right)$ is decided as follows. Values of physical quantity on sample points r₁, r₂, and r₃ around r₀ or P₁, P₂, and P₃ are approximated by Taylor expansions, which include derivatives $\partial P/\partial x, \partial P/\partial y,and\partial P/\partial z$ on r₀. These expansions can be considered as simultaneous equations, whose unknown quantities $\partial P/\partial x,\partial P/\partial y,and \partial P/\partial z$ can be solved by using r₁, r₂, and r₃ and P₁, P₂, and P₃. These expressions lead to matrix forms of $\left[\partial /\partial x\right],$$\left[\partial /\partial y\right], and \left[\partial /\partial \right]$. Higher order Taylor expansion and more P_is around r₀ are utilized.

A matrix form of $\left[{\nabla }^2\right]$is decided as follows. Poisson equation

  

$\begin{array}{c}{\nabla }^2\varphi =\rho \end{array}$(3)

is solved like

  

$\begin{array}{c}\varphi \left(r_1\right)=-\frac{1}{4\pi }\iiint \frac{\rho \left(r_2\right)}{r_{12}}d{v}_2.\end{array}$(4)

If we consider ϕ and ρ to be vectors $\overrightarrow{\varphi }$ and $\overrightarrow{\rho }$ like $\left[{\chi }_{ri}\right]$, (4) is expressed as a matrix operation

  

$\begin{array}{c}\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\\ \cdots \\ {\varphi }_n\end{array}\right]=-\frac{1}{4\pi }\left[\begin{array}{ccccc}\ddots & \frac{1}{r_{12}}& \frac{1}{r_{13}}& \cdots & \frac{1}{r_{1n}}\\ \frac{1}{r_{21}}& \ddots & \frac{1}{r_{23}}& \cdots & \frac{1}{r_{2n}}\\ \frac{1}{r_{31}}& \frac{1}{r_{32}}& \ddots & \cdots & \frac{1}{r_{3n}}\\ \cdots & \cdots & \cdots & \ddots & \cdots \\ \frac{1}{r_{n1}}& \frac{1}{r_{n2}}& \frac{1}{r_{n3}}& \cdots & \ddots \end{array}\right]\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\\ {\rho }_3\\ \cdots \\ {\rho }_n\end{array}\right]\end{array}$(5)

  

$or\overrightarrow{\varphi }=\left[R\right]\overrightarrow{\rho }$

where r_ij is the distance between r_i and r_j and diagonal elements of the matrix are decided analytically.

  

$\begin{array}{c}\left(5\right)canberewrittento{\left[R\right]}^{-1}\overrightarrow{\varphi }=\overrightarrow{\rho }. \end{array}$(6)

   Comparison (6) with (3) leads to

${\left[R\right]}^{-1}=\left[{\nabla }^2\right].$

Besides "well potential" like Figure 1 introduced in the preceding section is added to [F_ij] This electric potential has a negative value of W only inside a certain radius from the nucleus and adjusts the spreads of AOs.

Figure 1

 Typical well potential of depth W (in Hartree). The horizontal axis denotes R, the distance from the nucleus (in Bohr radius).

Thus, the Fock operator [F_ij] is found to be an entirely numerical matrix. Although [F_ij] isn't Hermitian at this stage, it becomes Hermitian if we change the element of the state vector from χ_ri to ${\chi }_{ri}\sqrt{{\omega }_i}$, where ω_i is the volume which r_i occupies. Thus, modified (2) can be solved by using a mathematical solver, even if A is an arbitrary vector potential. "HEEVR" of Intel Math Kernel Library^©, the double precision solver for complex Hermitian eigenvalue problems, is used in this study. In these calculations, 5774 sample points for the H atom and 8588 points for the 2nd row atoms were used. Calculations of corresponding matrices are heavy, but they can be done independently for all atoms constituting the molecule and total calculation time depends almost proportionally on the number of atoms.

Table 1 shows the AO energies of an argon atom, as examples of AO calculation by this method. They were compared with those calculated by the strict Hartree-Fock method [7]. AO #3-#5 and AO #7-#9 are degenerated, which suggests that they are corresponding to 2p and 3p orbitals, respectively. Figure 2 shows their shapes, which confirm these matchings.

Table 1

 AO energies of an argon atom, calculated by this method and strictly calculated counterparts [7] in Hartree.

Figure 2

 Shapes of AOs of argon atom. Blue and yellow designate plus and minus of AO, respectively.

Then obtained AOs were used to calculate molecular integrals and MOs were calculated by Roothaan's method.

In this method, the total energy of a molecule depends on the adjustable parameter α of Xα potential in (1). α for each atom in the molecule was adjusted to reproduce the total energy of the atom in its isolated situation. And number of sample points for each atom can affect the calculated total energy of the molecule. Table 2 displays the total energy of CO molecule with 3 kinds of sample point number per each atom, compared with the result by accurate Hartree Fock calculations [8] and demonstrates that adjustment of α and selection of sample point number were suitable.

Table 2

 Total energy of CO molecule in Hartree vs number of sample points used for each atom and the result by accurate Hartree Fock calculation [8].

3 Results

This method doesn't use any basis like GTO or STO, nevertheless calculated AOs of argon atom without A became complete spherical harmonic functions like Fig. 2. But if the potential loses spherical symmetry, for example by application of A, AOs fit to the potential and will become the best bases for the system. Although A is very weak generally, a small complex factor appears in AOs and it brings the magnetic character of molecules.

Figure 3

θ of 2s orbitals of carbon atom. The nucleus of an atom is located at (0,0), (10,0), (0,10), and (10,10) in atomic units. The value of θ is expressed by rainbow colors on sample points. The origin of A is located at (0,0).

A corresponding to uniform B₀ = (0,0,B₀) is not uniform but circulates origin, and magnetic quantities depend on the atom position. But the AO, which has complex factor at pointr (x, y, z)

  

$\begin{array}{c}{e}^{-\frac{i}{C}\theta \left(x,y,z\right)},where\theta \left(x,y,z\right)=B_0\frac{\left(xY-yX\right)}{2} \end{array}$(7)

and (X, Y, Z) is the position of the nucleus, forms physical quantities that aren't affected by the atom's position. This character is called "gauge invariance". (7) suggests that on the line connecting the origin and nucleus, θ ≡ 0. Locus of $\theta \left(x,y,z\right)\equiv const.$ is the line parallel to θ ≡ 0. The line for θ ˃ 0 lies at one side of θ ≡ 0 and the line for θ ˂ 0 at the other side.

Figure 3 shows θ of the 2s orbital of a carbon atom calculated by this method. The atom is located at 4 positions. These patterns indicate qualitatively that calculated θ follows (7). In the case of London-AO mentioned in the preceding section, complex factors are multiplied to real AO artificially. But in this method, suitable complex factor θ appears naturally in calculated AOs not only for uniform but also for arbitrary magnetic fields.

Uniform external magnetic field B₀ brings complex factors into MOs and generates current J as follows

  

$J=\frac{i}{2}\sum _k\left({\psi }_k\nabla {\psi }_k^{\text{*}}-{\psi }_k^{\text{*}}\nabla {\psi }_k\right).$(8)

J induces ∆B, which is easily calculated by Biot-Savart's law and ∆B modifies the magnetic field. The ratio $\sigma =-\left|∆B/B_0\right|$ at the position of the nucleus of an atom is called the "magnetic shielding constant" of that atom. Affluent data of σ is a useful measure to evaluate qualities of MOs related to magnetic character.

In a CH₃NH₂ molecule, σ(C), σ of carbon atom was calculated by this method as a function of W(C), the absolute value of its well potential depth, under the initial condition W(N) ≡ 0.400. As W has always a negative value, the absolute value of W is expressed just as W afterward. Results are shown in Figure 4. σ depends on W intensely and the maximum of σ is close to σ_exp, experimental σ. W_max(C) = 0.415, where σ(C) had its maximum, was decided by interpolation.

Figure 4

σ(C) vs absolute value of W(C) in CH₃NH₂ molecule. Horizontal bold line designates σ_exp.

Then σ(N) was calculated as a function of W(N) under the condition W(C) ≡ 0.415 decided at the preceding step. σ(N) had a maximum at W_max(N) = 0.342. Thus W(C) and W(N) were optimized alternatively and summarized in Table 3. W_max(N) of #4 almost becomes equal to that of #3 and the optimization was finished. Table 4 shows the final results under the #4 condition.

Table 3

 Optimized W of each step (in red) and W of partner atom optimized in the preceding step.

Figure 5 shows σs of various atoms in small molecules calculated in a similar way and corresponding σ_exps [9]. 5 bases for hydrogen atom and 9 bases for other atoms were used in these calculations.

Figure 5

 Magnetic shielding constants σs of atoms (molecule) calculated by this method and their experimental counterparts (hatched bars) [9]. In a C₃H₆ molecule, "C1" is the carbon atom to which 2 hydrogen atoms are bound, "C3" is the carbon atom to which 3 hydrogen atoms are bound, and "C2" is the central carbon atom.

Table 4

Finally calculated σ(C) and σ(N) of CH₃NH₂ molecule and their experimental counterparts [9].

Here σs of 3 carbon atoms of propylene (C₃H₆) were obtained by optimizing W of each atom independently because 3 atoms have different bonding with their adjacent atoms. This treatment improves coincidence with experimental counterparts. But as to the same atoms in the same circumstance, the same AOs can be used. The obtained results in Table 4 and Figure 5 agreed fairly good with their experimental counterparts.

If the gauge invariance of MO is complete, calculated σs will not depend on the position of the molecule at all. But its incompleteness remains considerable. σs of atoms in a CH₃NH₂ molecule calculated in conditions that it was moved parallelly or rotated were summarized in Table 5. Atoms are numbered as in Figure 6.

Figure 6

 Atom numberings in CH₃NH₂ molecule.

Table 5

"b" in the left-most column means that σs in the same row were calculated by the method of this study in the condition that the carbon atom was located at the origin of A. "a" or "c" means that the molecule was moved by 3.0Å to minus x-direction or plus x direction respectively. "x" means that the molecule was rotated clockwise around the x-axis by 45.0 degrees.

On the other hand, Table 6 shows σs calculated by using AOs by solving the HFS equation which doesn't include magnetic interaction terms. Here σs vary extremely by moving the molecule and the effect of improved AOs is demonstrated. To minimize errors caused by incomplete gauge invariance, it was very effective to locate the centers of the molecules at the origin of A. Calculations in Table 4 and Figure 5 were done by this way. But the effect of this method is limited for small molecules and improvement of the accuracy of gauge invariance is strongly required.

Table 6

σs calculated by using AOs without considering magnetic interaction as bases. "a", "b", "c" and "x" have the same meaning as Table 5.

4 Discussion

The coincidence of the maximum of σ and σ_exp is very useful for obtaining accurate σ. The author speculates its mechanism as follows. B₀ generates J and J has self-interaction energy E_J, which is the sum of the contribution from each atom (designated as A) like,

  

$\begin{array}{c}{E}_J=\sum _A{E}_J^{\text{A}}=\sum _A\iiint {}^{{}^{{}^A}}\Delta B\cdot \left(r\times J\right)dv\end{array}$(9)

where (r × J) is the magnetic moment of J and ∆B is the magnetic field that J generated. ∆B is a function of r but can be approximated by −σ(A)B₀ in the area belonging to atom A and $E_J^A$ is represented like

  

$\begin{array}{c}{E}_J^{\text{A}}\approx -\sigma \left(\text{A}\right)\underset{\_}{B_0\cdot \iiint {}^{{}^{{}^A}}\left(r\times J\right)dv}. \end{array}$(10)

When W(A) is increased gradually, MO around atom A will have its optimum shape locally and $E_J^{\text{A}}$ will have its minimum at some W(A), W_opt(A). If the underlined part of (10) changes far more slowly than σ(A) and is positive, σ(A) will have its maximum and come near σ_exp(A) at W_opt(A). In the case that the underlined part is negative, σ(A) will have its minimum. Such cases sometimes occured in optimization steps.

$E_J^{\text{A}}$s in (9) are not completely independent of each other and W of one atom may affect other atoms to some extent. Therefore SCF-like optimization was necessary to decide optimized Ws.

In the case of other physical quantities like molecular magnetism, situation will be different from σ. As to σ, (10) holds and the maximum of σ comes near σ_exp at the optimized W. But this lucky situation will not be expected for other physical quantities.

5 Conclusion and Future Challenges

A new method to calculate MOs of molecules in an arbitrary magnetic field by the LCAO-MO method has been developed. Calculated magnetic shielding constants demonstrated that this method works well at least in the case of uniform magnetic fields and small molecules and clarified that careful preparation of bases of the LCAO-MO method is very effective to obtain better physical quantity. The author is searching now new applications where a nonuniform magnetic field is essential and/or carefully tuned-up bases are required.

On the other hand, problems impeding a wider range of applications of this method remain. Accuracy of gauge invariance of AOs is insufficient. Calculation of numerical eigenvalue problems is too time-consuming because of the large number of sample points and the poor parallel efficiency of its solver.

References

• [1]   F. London, J. Phys. Radium, 8, 397 (1937). doi:10.1051/jphysrad:01937008010039700
• [2]   K. Nakagawa, J. Comput. Chem. Jpn., 14, 50 (2015). doi:10.2477/jccj.2015-0025
• [3]   J. C. Slater, Quantum Theory of Molecules and Solids Vol.4 Magraw-Hill(1974)
• [4]   D. E. Ellis, H. Adachi, F. W. Averill, Surf. Sci., 58, 497 (1976). doi:10.1016/0039-6028(76)90485-4
• [5]   F. W. Averill, D. E. Ellis, J. Chem. Phys., 59, 6412 (1973). doi:10.1063/1.1680020
• [6]   https://handwiki.org/wiki/Physics:Atomic_units
• [7]   S. Fraga, et al., Handbook of Atomic Data, Elsevier (1976)
• [8]   P. A. Christiansen, E. A. McCullough, Jr., J. Chem. Phys., 67, 1877 (1977). doi:10.1063/1.435144
• [9]   A. M. Teale, O. B. Lutnæs, T. Helgaker, D. J. Tozer, J. Gauss, J. Chem. Phys., 138, 024111 (2013). doi:10.1063/1.4773016 PMID:23556715