LoProp charges and polarisabilities

LoProp charges and polarisabilities#

The localised properties (LoProp) method [GLKarlstrom04] is a computationally inexpensive scheme that gives physically meaningful localised charges and polarisabilities that sum up to the corresponding molecular properties.

In this scheme, we classify atomic orbitals in the basis set to be either occupied or virtual based on the associated atomic electron configurations. The former category includes both fully and partially occupied atomic sub-shells.

We note that the family of atomic natural orbital (ANO) basis sets is well suited for LoProp calculations as the basis functions are close to the orbitals of the atoms.

Reference LoProp calculations#

As an example, let us consider hydrogen fluoride at the Hartree–Fock/6-31G level of theory. We first run a reference LoProp calculation with the built-in module in the VeloxChem program and will thereafter repeat the LoProp calculation step-by-step in the notebook.

import veloxchem as vlx
import numpy as np
from numpy.linalg import inv

np.set_printoptions(precision=4, suppress=True, linewidth=132)

mol_str = """
F        0.00000000    0.00000000     -0.1733
H        0.00000000    0.00000000      1.5597
"""

molecule = vlx.Molecule.read_molecule_string(mol_str, units='au')
basis = vlx.MolecularBasis.read(molecule, "6-31G")

natoms = molecule.number_of_atoms()
norb = basis.get_dimensions_of_basis(molecule)

# SCF optimization
scf_drv = vlx.ScfRestrictedDriver()
scf_results = scf_drv.compute(molecule, basis)

# LoProp
loprop_drv = vlx.LoPropDriver()
loprop_out = loprop_drv.compute(molecule, basis, scf_results)

print("Localized charges (a.u.):")
print(f"F: {loprop_out['localized_charges'][0] : .4f}")
print(f"H: {loprop_out['localized_charges'][1] : .4f}")

print("\nLocalized polarizabilities (a.u.):")
print("     xx      xy      xz      yy      yz      zz")
print(f"F: {loprop_out['localized_polarizabilities'][0]}")
print(f"H: {loprop_out['localized_polarizabilities'][1]}")

Localized charges (a.u.):
F: -0.1191
H:  0.1191

Localized polarizabilities (a.u.):
     xx      xy      xz      yy      yz      zz
F: [ 0.4371  0.     -0.      0.4371 -0.      1.917 ]
H: [ 0.226  -0.     -0.      0.226  -0.      1.9525]

LoProp transformation#

A LoProp transformation is performed of the atomic orbitals to obtain an orthonormal set of basis functions

\[ | \boldsymbol{\chi}^\mathrm{L} \rangle = | \boldsymbol{\chi} \rangle \mathbf{T} \]

This transformation is decomposed into five steps

\[ \mathbf{T} = \mathbf{T}_0 \times \mathbf{T}_1 \times \mathbf{T}_2 \times \mathbf{T}_3 \times \mathbf{T}_4 \]

T0: Reordering of atomic orbitals
T1: Gram–Schmidt orthonormalisation for each atom block
T2: Lödwin orthonomalisation on occupied and virtual
T3: occupied to virtual projection
T4: Lödwin virtual

with the overall aim to retain the atomic character of the occupied orbitals.

We start with the overlap matrix, \(S\), in the original AO basis:

S = scf_results['S']

print(S)

[[1.     0.2443 0.1716 0.0341 0.0601 0.     0.     0.     0.     0.     0.    ]
 [0.2443 1.     0.7635 0.2275 0.3321 0.     0.     0.     0.     0.     0.    ]
 [0.1716 0.7635 1.     0.4306 0.6374 0.     0.     0.     0.     0.     0.    ]
 [0.0341 0.2275 0.4306 1.     0.6583 0.     0.     0.2409 0.5971 0.     0.    ]
 [0.0601 0.3321 0.6374 0.6583 1.     0.     0.     0.1098 0.4105 0.     0.    ]
 [0.     0.     0.     0.     0.     1.     0.4995 0.     0.     0.     0.    ]
 [0.     0.     0.     0.     0.     0.4995 1.     0.     0.     0.     0.    ]
 [0.     0.     0.     0.2409 0.1098 0.     0.     1.     0.4995 0.     0.    ]
 [0.     0.     0.     0.5971 0.4105 0.     0.     0.4995 1.     0.     0.    ]
 [0.     0.     0.     0.     0.     0.     0.     0.     0.     1.     0.4995]
 [0.     0.     0.     0.     0.     0.     0.     0.     0.     0.4995 1.    ]]

T0: Atomic orbital rearrangement#

../../_images/loprop_T0.svg — Fig. 10 Re-arranging atomic orbitals.#

For hydrogen floride with the user defined order of atoms (F before H), the AOs are arranged as follow:

\(\chi_{1s}^{F}\), \(\chi_{2s}^{F}\), \(\chi_{3s}^{F}\) \(\chi_{1s}^{H}\), \(\chi_{2s}^{H}\) \(\chi_{1py}^{F}\), \(\chi_{2py}^{F}\), \(\chi_{1pz}^{F}\), \(\chi_{2pz}^{F}\), \(\chi_{1px}^{F}\), \(\chi_{2px}^{F}\)

where orbitals classified as occupied are high-lighted in red. The \(T_0\) transformation re-arranges the orbitals as follows:

\(\chi_{1s}^{F}\), \(\chi_{2s}^{F}\), \(\chi_{3s}^{F}\), \(\chi_{1px}^{F}\), \(\chi_{1py}^{F}\),\(\chi_{1pz}^{F}\), \(\chi_{2px}^{F}\), \(\chi_{2py}^{F}\),\(\chi_{2pz}^{F}\), \(\chi_{1s}^{H}\), \(\chi_{2s}^{H}\)

rearranged_indices=[]

for atom in range(natoms):
    indices, angmoms = vlx.get_basis_function_indices_for_atom(molecule, basis, atom)
    rearranged_indices.append(indices)
    
rearranged_indices = np.concatenate(rearranged_indices)

print('Rearranged AO indices：', rearranged_indices)

Rearranged AO indices： [ 0  1  2  9  5  7 10  6  8  3  4]

T0 = np.zeros((norb, norb))

for col, row in enumerate(rearranged_indices):
    T0[row, col] = 1

S0 = np.einsum('ba,bc,cd->ad', T0, S, T0)
print(S0)

[[1.     0.2443 0.1716 0.     0.     0.     0.     0.     0.     0.0341 0.0601]
 [0.2443 1.     0.7635 0.     0.     0.     0.     0.     0.     0.2275 0.3321]
 [0.1716 0.7635 1.     0.     0.     0.     0.     0.     0.     0.4306 0.6374]
 [0.     0.     0.     1.     0.     0.     0.4995 0.     0.     0.     0.    ]
 [0.     0.     0.     0.     1.     0.     0.     0.4995 0.     0.     0.    ]
 [0.     0.     0.     0.     0.     1.     0.     0.     0.4995 0.2409 0.1098]
 [0.     0.     0.     0.4995 0.     0.     1.     0.     0.     0.     0.    ]
 [0.     0.     0.     0.     0.4995 0.     0.     1.     0.     0.     0.    ]
 [0.     0.     0.     0.     0.     0.4995 0.     0.     1.     0.5971 0.4105]
 [0.0341 0.2275 0.4306 0.     0.     0.2409 0.     0.     0.5971 1.     0.6583]
 [0.0601 0.3321 0.6374 0.     0.     0.1098 0.     0.     0.4105 0.6583 1.    ]]

T1: Gram–Schmidt orthogonalization of atomic blocks#

../../_images/loprop_T1.svg — Fig. 11 Diagonalization of the atomic blocks.#

Orthonormality in the atomic blocks of the overlap matrix is reached with a Gram–Schmidt orthogonalization, or equivalently with a Cholesky factorization. Using Gram-Schmidt orthogonalization ensures that the atomic subspace of the occupied basis functions is left intact. We have

\(S_0 = L L^T\)

where \(L\) is a lower triangular matrix. We get

\(I = L^{-1} S_0 [L^T]^{-1}\)

and identify

\(T_1 = [L^T]^{-1}\)

T1 = np.zeros((norb,norb))

norb_per_atom, indices_occ, indices_virt = loprop_drv.get_ao_indices(molecule, basis)

print('Occupied orbitals:', indices_occ)
print('Virtual orbitals:', indices_virt, '\n')

ri = 0
for atom in range(natoms):

    rf = ri + norb_per_atom[atom]

    L = np.linalg.cholesky(S0[ri:rf, ri:rf])
    T1[ri:rf, ri:rf] = np.linalg.inv(L.T)

    ri += norb_per_atom[atom]

S1 = np.einsum('ba,bc,cd->ad', T1, S0, T1)
print(S1)

Occupied orbitals: [0, 1, 3, 4, 5, 9]
Virtual orbitals: [2, 6, 7, 8, 10] 

[[ 1.     -0.      0.      0.      0.      0.      0.      0.      0.      0.0341  0.0501]
 [ 0.      1.     -0.      0.      0.      0.      0.      0.      0.      0.226   0.2372]
 [ 0.     -0.      1.      0.      0.      0.      0.      0.      0.      0.3975  0.4416]
 [ 0.      0.      0.      1.      0.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      1.      0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      1.      0.      0.      0.      0.2409 -0.0647]
 [ 0.      0.      0.      0.      0.      0.      1.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      1.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      1.      0.5504  0.0641]
 [ 0.0341  0.226   0.3975  0.      0.      0.2409  0.      0.      0.5504  1.     -0.    ]
 [ 0.0501  0.2372  0.4416  0.      0.     -0.0647  0.      0.      0.0641 -0.      1.    ]]

T2: Lödwin orthonormalization of occupied and virtual blocks#

../../_images/loprop_T2.svg — Fig. 12 Diagonalization of the occupied and virtual blocks.#

The use of Löwdin orthonormalization ensures that the orbitals in the occupied and virtual subspaces are made orthonormal with a minimal modification as compared to the original atomic orbitals.

T2 = np.zeros((norb, norb))

S1_oo = S1[indices_occ, :][:, indices_occ]
T2_oo = loprop_drv.lowdin_orthonormalize(S1_oo)

nocc = len(indices_occ)
nvirt = len(indices_virt)

for r in range(nocc):
    for c in range(nocc):
        T2[indices_occ[r], indices_occ[c]] = T2_oo[r, c]

S1_vv = S1[indices_virt, :][:, indices_virt]
T2_vv = loprop_drv.lowdin_orthonormalize(S1_vv)

for r in range(nvirt):
    for c in range(nvirt):
        T2[indices_virt[r], indices_virt[c]] = T2_vv[r, c]
        
S2 = np.einsum('ba,bc,cd->ad', T2, S1, T2)
print(S2)

[[ 1.     -0.     -0.0209 -0.      0.     -0.      0.      0.     -0.0121  0.      0.0573]
 [-0.      1.     -0.1143 -0.      0.     -0.      0.      0.     -0.0765 -0.      0.2767]
 [-0.0209 -0.1143  1.     -0.      0.     -0.0414  0.      0.      0.      0.4636 -0.    ]
 [-0.     -0.     -0.      1.     -0.     -0.      0.      0.     -0.     -0.     -0.    ]
 [ 0.      0.      0.     -0.      1.      0.      0.      0.      0.      0.     -0.    ]
 [-0.     -0.     -0.0414 -0.      0.      1.      0.      0.     -0.0698 -0.     -0.0505]
 [ 0.      0.      0.      0.      0.      0.      1.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      1.      0.      0.      0.    ]
 [-0.0121 -0.0765  0.     -0.      0.     -0.0698  0.      0.      1.      0.5822 -0.    ]
 [ 0.     -0.      0.4636 -0.      0.     -0.      0.      0.      0.5822  1.     -0.1497]
 [ 0.0573  0.2767 -0.     -0.     -0.     -0.0505  0.      0.     -0.     -0.1497  1.    ]]

T3: Occupied–virtual projection#

../../_images/loprop_T3.svg — Fig. 13 Projecting out the occupied space from the virtual.#

Components from the occipied subspace are projected out of the virtual subspace by means of Gram–Schmidt orthonormalizaton.

T3 = np.identity(norb)

S2_ov = S2[indices_occ, :][:, indices_virt]

for i, index_occ in enumerate(indices_occ):
    for j, index_virt in enumerate(indices_virt):
        T3[index_occ, index_virt] = - S2_ov[i, j]

S3 = np.einsum('ba,bc,cd->ad', T3, S2, T3)
print(S3)

[[ 1.     -0.     -0.     -0.      0.     -0.      0.      0.     -0.      0.      0.    ]
 [-0.      1.      0.     -0.      0.     -0.      0.      0.      0.     -0.      0.    ]
 [-0.      0.      0.7698  0.     -0.     -0.      0.      0.     -0.2818  0.      0.1001]
 [-0.     -0.      0.      1.     -0.     -0.      0.      0.      0.     -0.      0.    ]
 [ 0.      0.     -0.     -0.      1.      0.      0.      0.     -0.      0.      0.    ]
 [-0.     -0.     -0.     -0.      0.      1.      0.      0.     -0.     -0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      1.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.      0.      1.      0.      0.      0.    ]
 [-0.      0.     -0.2818  0.     -0.     -0.      0.      0.      0.6502  0.      0.1055]
 [ 0.     -0.      0.     -0.      0.     -0.      0.      0.      0.      1.     -0.    ]
 [ 0.      0.      0.1001  0.      0.      0.      0.      0.      0.1055 -0.      0.8952]]

T4: Lödwin orthogonalizatio of virtual blocks#

../../_images/loprop_T4.svg — Fig. 14 Diagonalization the virtual subspace.#

The virtual subspace is Löwdin orthonmalised.

T4 = np.identity(norb)

S3_vv = S3[indices_virt, :][:, indices_virt]
T4_vv = loprop_drv.lowdin_orthonormalize(S3_vv)

for i, index_virt1 in enumerate(indices_virt):
    for j, index_virt2 in enumerate(indices_virt):
        T4[index_virt1, index_virt2] = T4_vv[i, j]
        
S4 = np.einsum('ba,bc,cd->ad', T4, S3, T4)
print(S4)

[[ 1. -0. -0. -0.  0. -0.  0. -0. -0.  0.  0.]
 [-0.  1.  0. -0.  0. -0.  0.  0.  0. -0.  0.]
 [-0.  0.  1.  0. -0. -0.  0. -0.  0.  0. -0.]
 [-0. -0.  0.  1. -0. -0.  0. -0.  0. -0.  0.]
 [ 0.  0. -0. -0.  1.  0. -0. -0. -0.  0.  0.]
 [-0. -0. -0. -0.  0.  1.  0. -0. -0. -0.  0.]
 [ 0.  0.  0.  0. -0.  0.  1.  0.  0.  0.  0.]
 [-0.  0. -0. -0. -0. -0.  0.  1.  0.  0. -0.]
 [-0.  0.  0.  0. -0. -0.  0.  0.  1.  0. -0.]
 [ 0. -0.  0. -0.  0. -0.  0.  0.  0.  1. -0.]
 [ 0.  0. -0.  0.  0.  0.  0. -0. -0. -0.  1.]]

The complete LoProp transformation matrix is formed.

T = np.einsum('ab,bc,cd,de,ef->af', T0, T1, T2, T3, T4)

The orbitals in the LoProp basis are atomic in character:

print('   F 1s    F 2s    F 2px   F 2py   F 2pz   H 1s\n', '-'*49)
print(T[:,indices_occ])

   F 1s    F 2s    F 2px   F 2py   F 2pz   H 1s
 -------------------------------------------------
[[ 0.9997 -0.2541  0.      0.     -0.0022  0.0123]
 [ 0.0032  1.0527 -0.      0.      0.0229 -0.1253]
 [ 0.      0.      0.      0.      0.      0.    ]
 [-0.0183 -0.1215 -0.      0.     -0.1295  1.045 ]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.     -0.      1.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.0033  0.0222 -0.      0.      1.0237 -0.1295]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.     -0.      1.     -0.     -0.     -0.    ]
 [ 0.      0.      0.      0.      0.      0.    ]]

Localization of first-order properties#

A molecular first-order property associated with an expctation value of a one-electron operator can be expressed as

\[ \langle \hat{\Omega} \rangle = \mathrm{tr} \big[ \Omega D \big] \]

where \(\Omega\) and \(D\) are the property and one-particle reduced density matrices in original AO basis, respectively. This trace can equally well be evaluated in the LoProp basis after transforming the matrices.

From the LCAO expansions of MOs

\[ | \boldsymbol{\phi} \rangle = | \boldsymbol{\chi} \rangle C = | \boldsymbol{\chi}^\mathrm{L} \rangle T^{-1} C \]

we identify

\[ C^\mathrm{L} = T^{-1} C \]

and consequently

\[ D^\mathrm{L} = T^{-1} D \, [T^\dagger]^{-1} = T^{-1} D \, [T^{-1}]^\dagger \]

The transformation of property matrices is straightforwardly shown to be given by

\[ \Omega^\mathrm{L} = T^\dagger \Omega T \]

We arrive at an decomposition of first-order properties in atomic contributions according to

\[ \langle \hat{\Omega} \rangle = \sum_{A} \Omega^\mathrm{L}_{A}; \qquad \Omega^\mathrm{L}_{A} = \sum_{B} \Omega^\mathrm{L}_{AB}; \qquad \Omega^\mathrm{L}_{AB} = \sum_{\alpha \in A} \sum_{\beta \in B} \Omega^\mathrm{L}_{\alpha\beta} D^\mathrm{L}_{\beta\alpha} \]

where \(A\) and \(B\) are atom indices and \(\alpha\) and \(\beta\) are atomic orbital indices.

def atomic_property(Omega, D):
    
    Omega_AB = np.zeros((natoms, natoms))

    Omega_L = np.einsum('ba,bc,cd->ad', T, Omega, T)
    D_L = np.einsum('ab,bc,dc->ad', inv(T), D, inv(T))
    
    idx_A = 0
    for A in range(natoms):

        idx_B = 0
        for B in range(natoms):
            
            P = Omega_L[idx_A : idx_A + norb_per_atom[A], idx_B : idx_B + norb_per_atom[B]]
            Q = D_L[idx_A : idx_A + norb_per_atom[A], idx_B : idx_B + norb_per_atom[B]]
        
            Omega_AB[A,B] = np.matmul(P.T, Q).trace()
                    
            idx_B += norb_per_atom[B]
            
        idx_A += norb_per_atom[A]
        
    return Omega_AB

Localized charges#

The total electronic charge is equal to

\[ Q = -e \, \mathrm{tr} \big[ S D \big] \]

The atomic charges are obtained from the partial summations over atomic orbitals in the orthonormal LoProp basis

\[ Q_A = - e \sum_{\alpha \in A} D^\mathrm{L}_{\alpha\alpha} \]

# atomic numbers and therefore nuclear charge in atomic units
Z_A = molecule.elem_ids_to_numpy()
D = 2 * scf_results['D_alpha']

Q_A = - atomic_property(S, D).sum(axis = 1)
Q = Q_A.sum()

print(f"Total electronic charge: {Q : .4f}")
print(f"Localized atomic charges: {Q_A + Z_A}")

Total electronic charge: -10.0000
Localized atomic charges: [-0.1191  0.1191]

Localized dipole moments#

With \(\hat{\Omega} = -e \, \hat{z}\), the associated property is the electronic dipole moment along the bond axis.

The total molecular electronic dipole moment is equal to

\[ \langle \hat{\mu} \rangle = \mathrm{tr} \big[ \mu D \big] \]

dipole_drv = vlx.ElectricDipoleIntegralsDriver()
dipole_mats = dipole_drv.compute(molecule, basis)

mu = - dipole_mats.z_to_numpy()

mu_AB = atomic_property(mu, D)

print(f'Total dipole moment: {mu_AB.sum() : .4f}')
print('\nmu_AB:\n', mu_AB)

Total dipole moment:  0.9033

mu_AB:
 [[ 2.0947  0.1471]
 [ 0.1471 -1.4855]]

We note that the nuclear contribution to the total molecular dipole moment here vanished as due to our choice of placing the origin of the coordinate system at the center of nuclear charge.

The local dipole moment is origin dependent and local origins are therefore used.

../../_images/loprop_local_origin.svg — Fig. 15 Bond midpoint coordinates.#

We introduce

\[ \boldsymbol{R}_{AB} = \frac{\boldsymbol{R}_{A} + \boldsymbol{R}_{B}}{2} \]

as the bond midpoint coordinate. When \(\boldsymbol{R}_{A} = \boldsymbol{R}_{B}\), it reduces to the atomic coordinate.

R = np.zeros((natoms, natoms, 3))

coords = molecule.get_coordinates()

for A in range(natoms):
    for B in range(natoms):
    
        R[A,B,:] = (coords[A,:] + coords[B,:]) / 2

With use of local origins, we define a gauge-orgin independent dipole moment

\[ \boldsymbol{\mu}^\mathrm{L}_{(AB)} = -e \sum_{\alpha \in A} \sum_{\beta \in B} \langle \chi^\mathrm{L}_\alpha | (\hat{\boldsymbol{r}} - \boldsymbol{R}_{AB}) | \chi^\mathrm{L}_\beta \rangle D^\mathrm{L}_{\beta\alpha} \]

It can be determined from the relation

\[ \boldsymbol{\mu}^\mathrm{L}_{(AB)} = \boldsymbol{\mu}^\mathrm{L}_{AB} - Q_A \boldsymbol{R}_{A} \delta_{AB} \]

print('mu_(AB):\n', mu_AB - np.diag(Q_A) * R[:,:,2])

mu_(AB):
 [[ 0.5144  0.1471]
 [ 0.1471 -0.1116]]

Localization of second-order properties#

Polarizabilities from response functions#

The components of the electric-dipole polarizability tensor ia calculated from the linear response function

\[ \alpha(-\omega;\omega) = - \langle\langle \hat{\mu}; \hat{\mu} \rangle\rangle_\omega = -{\mu^{[1]}}^\dagger N(\omega) \]

where we have omitted the Cartesian tensor indices associated with the components of the dipole operators. For illustration, we will focus on the \(zz\)-component in the static limit, \(\omega = 0\).

The strucure of the property gradient is

\[\begin{split} \mu^{[1]} = \begin{pmatrix} g \\ -g^* \\ \end{pmatrix} \end{split}\]

where the elements of \(g\) refer to pairs of occupied–virtual orbital indices. With use of normalized, spin-adapted, electron excitation operators

\[ \hat{E}_{si}^\dagger = \frac{1}{\sqrt{2}} \big( \hat{a}^\dagger_{s, \alpha} \hat{a}_{i, \alpha} + \hat{a}^\dagger_{s, \beta} \hat{a}_{i, \beta} \big) \]

we get

\[ g_{si} = \langle 0 | \big[ \hat{E}_{si}, \hat{\mu} \big] | 0 \rangle = \sqrt{2} \mu_{si} \]

The response vector \(N(\omega)\) is introduced as the solution vector to the linear response equation

\[ N(\omega) = \Big( E^{[2]} - \hbar \omega S^{[2]} \Big)^{-1} \mu^{[1]} \]

The structure of the response vectors is

\[\begin{split} N(\omega) = \begin{pmatrix} Z \\ -Y^* \\ \end{pmatrix} \end{split}\]

and in the static limit, we have \(Y = Z\).

lrs_drv = vlx.LinearResponseSolver()
lrs_out = lrs_drv.compute(molecule, basis, scf_results)

alpha = - lrs_out['response_functions'][('z', 'z', 0)]

print(f"Molecular polarizability (zz-component): {alpha : .4f} a.u.")

Molecular polarizability (zz-component):  3.8695 a.u.

Polarizabilities from perturbed density matrices#

The upper and lower parts of the response vector can be unpacked into the occupied–virtual (ov) and virtual–occupied (vo) blocks of associate MO matrices according to

\[\begin{split} \kappa = \begin{pmatrix} 0 & Z \\ Y^* & 0 \\ \end{pmatrix} \end{split}\]

def vec2mat(N, nocc, norb):
    
    n = nocc * (norb - nocc)
    
    kappa = np.zeros((norb,norb))
    
    kappa[:nocc, nocc:norb] = N[:n].reshape(nocc, norb - nocc)
    kappa[nocc:, :nocc] = - N[n:].reshape(nocc, norb - nocc).T
    
    return kappa

The polarizability can then be evaluated in the orignal AO basis according to

\[ \alpha(-\omega; \omega) = \mathrm{tr}\big[ \mu \,\kappa \big] = \sum_{\alpha,\beta} \mu_{\alpha\beta} D^{(1)}_{\alpha\beta} \]

where we have introduced the perturbed density matrix

\[ D^{(1)}_{\alpha\beta} = \sum_{p,q} C^*_{\alpha p} C_{\beta q} \kappa_{qp} \]

It can also be exressed in the LoProp basis

\[ \alpha(-\omega; \omega) = \sum_{A,B} \alpha^\mathrm{L}_{AB}; \qquad % \alpha^\mathrm{L}_{AB} = \sum_{\alpha \in A} \sum_{\beta \in B} \mu^\mathrm{L}_{\alpha\beta} D^{\mathrm{L},(1)}_{\beta\alpha} \]

nocc = molecule.number_of_alpha_electrons()

solution_z_0 = lrs_drv.get_full_solution_vector(lrs_out['solutions'][('z', 0)])

N = - solution_z_0

kappa = vec2mat(N, nocc, norb)

C = scf_results['C_alpha']
# factor of sqrt(2) due to use of normalized spin-adapted excitation operators
# in the formation of the property gradient, mu[1].
D1 = np.einsum('QP,aP,bQ->ab', kappa, C, C) * np.sqrt(2)

alpha_AB = atomic_property(mu, D1)

print(f"alpha = {alpha_AB.sum() : .6f} a.u.")
print('\nalpha_AB:\n', alpha_AB)

alpha =  3.869484 a.u.

alpha_AB:
 [[0.6304 0.2307]
 [0.2307 2.7776]]

Localized polarizabilities#

We define gauge-origin independent local polarizability

\[ \alpha^\mathrm{L}_{(AB)} = -e \sum_{\alpha \in A} \sum_{\beta \in B} \langle \chi^\mathrm{L}_\alpha | (\hat{\boldsymbol{r}} - \boldsymbol{R}_{AB}) | \chi^\mathrm{L}_\beta \rangle D^{\mathrm{L},(1)}_{\beta\alpha} \]

It can be determined from the relation

\[ \alpha^\mathrm{L}_{(AB)} = \alpha^\mathrm{L}_{AB} - \Delta Q_A \boldsymbol{R}_{A} \delta_{AB} \]

where the external field derivative of the change in charge on atom \(A\) is

\[ \Delta Q_A = -e \sum_{\alpha \in A} D^{\mathrm{L},(1)}_{\alpha\alpha} \]

As charge is preserved, we have

\[ \sum_A \Delta Q_A = 0 \]

dQ_A = - atomic_property(S, D1).sum(axis = 1)

# Lower case atomic indices for gauge-origin independent properties
alpha_ab = alpha_AB - np.diag(dQ_A) * R[:,:,2]

print('dQ_A:\n', dQ_A)
print('alpha_(AB):\n', alpha_ab)

dQ_A:
 [-1.5231  1.5231]
alpha_(AB):
 [[0.3664 0.2307]
 [0.2307 0.4019]]

Next, we define an auxiliary charge transfer matrix with elements, \(\Delta Q_{AB}\), as the external field derivative of the charge transfer from center \(A\) to \(B\). Although this quantity is not uniquely defined, it has to fulfill

\[ \sum_{A \neq B} \Delta Q_{AB} = \Delta Q_{A} \]

and

\[ \Delta Q_{BA} = - \Delta Q_{AB} \]

In the LoProp scheme, it assumed that \(\Delta Q_{AB}\) should be as small as possible and also that the charge transfer is short range. To accommodate these desired properties, a Lagrangian is formed

\[ \mathcal{L}(\Delta Q_{AB}, \lambda_A) = \sum_{AB} \big[ \Delta Q_{AB} \big]^2 f(\Delta R_{AB}) + \sum_{A} \lambda_A \big( \sum_{A \neq B} \Delta Q_{AB} - \Delta Q_A \big) \]

where \(\Delta R_{AB} = |\boldsymbol{R}_A - \boldsymbol{R}_B|\) and \(f\) is a penatly function preventing long-distance charge transfer chosen as

\[ f(r) = \exp\Big(\alpha \Big(\frac{r}{R_A^\mathrm{BS} + R_B^\mathrm{BS}}\Big)^2\Big) \]

The constant \(\alpha\) is chosen equal to \(2.0\) and \(R_A^\mathrm{BS}\) is the Bragg–Slater radius of atom \(A\).

def f(A, B):
    
    # Bragg-Slater radii for F and H
    bohr2ang = 0.529177249
    R_BS = np.array([0.50, 0.25]) / bohr2ang
    
    # interatomic distance
    dR_AB = np.linalg.norm(R[A,A,:] - R[B,B,:])
        
    return np.exp(2.0 * (dR_AB / (R_BS[A] + R_BS[B]))**2)

The stationary point of the Lagrangian allows for a determination of the Lagrange multipliers from the equation

\[ L \lambda = \Delta Q \]

where

\[ L_{AB} = \frac{1}{2 f(\Delta R_{AB})} + C \]

and

\[ L_{AA} = -\sum_{A \neq B} \frac{1}{2 f(\Delta R_{AB})} + C \]

The constant \(C\) is determined from the requirement of charge conservation, or zero net charge transfer.

Once the Lagrange multipliers are determined, we have

\[ \Delta Q_{AB} = - \frac{\lambda_{A} - \lambda_{B}}{2 f(\Delta R_{AB})} \]

L_AB = np.zeros((natoms, natoms))

for A in range(natoms):
    for B in range(natoms):
        L_AB[A,B] = 0.5 / f(A,B)

    L_AB[A,A] -= sum(L_AB[A,:])

L_AB = L_AB + 2.0 * np.max(np.abs(L_AB))

lambda_A = np.linalg.solve(L_AB, dQ_A)

print(f'lambda_A: {lambda_A}')

lambda_A: [ 30.2963 -30.2963]

dQ_AB = np.zeros((natoms,natoms))

for A in range(natoms):
    for B in range(natoms):

        dQ_AB[A,B] = - 0.5 * (lambda_A[A] - lambda_A[B]) / f(A,B)

print(dQ_AB)

[[-0.     -1.5231]
 [ 1.5231 -0.    ]]

The \(zz\)-component of the bond polarizability is defined by

\[ \alpha^\mathrm{bond}_{AB} = \Delta Q_{AB} \Delta \boldsymbol{R}_{AB, z} \]

alpha_bond_AB = np.zeros((natoms, natoms))

for A in range(natoms):
    for B in range(natoms):
        alpha_bond_AB[A,B] = dQ_AB[A,B] * (R[A,A,2] - R[B,B,2])
        
print(alpha_bond_AB)

[[-0.      2.6396]
 [ 2.6396 -0.    ]]

The LoProp polarizability for atom \(A\) is obtained as

\[ \alpha^\mathrm{L}_A = \sum_B \Big( \alpha^\mathrm{L}_{(AB)} + \frac{1}{2} \alpha^\mathrm{bond}_{AB} \Big) \]

and it fulfills

\[ \sum_A \alpha^\mathrm{L}_A = \alpha(-\omega;\omega) \]

alpha_A = (alpha_ab + 0.5 * alpha_bond_AB).sum(axis = 1)

print(f"LoProp atomic polarizabilities: {alpha_A}")
print(f"Molecular polarizability: {alpha_A.sum() : .4f} a.u.")

LoProp atomic polarizabilities: [1.917  1.9525]
Molecular polarizability:  3.8695 a.u.