Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

LoProp charges and polarisabilities

The localised properties (LoProp) method Gagliardi et al. (2004) 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/ANO-S-VDZP 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=200)

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, "ANO-S-VDZP")

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

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

# LoProp
loprop_drv = vlx.PEForceFieldGenerator()
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

∣χL⟩=∣χ⟩T| \boldsymbol{\chi}^\mathrm{L} \rangle = | \boldsymbol{\chi} \rangle \mathbf{T}
(1)

This transformation is decomposed into five steps

T=T0×T1×T2×T3×T4\mathbf{T} = \mathbf{T}_0 \times \mathbf{T}_1 \times \mathbf{T}_2 \times \mathbf{T}_3 \times \mathbf{T}_4
(2)

  • 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, SS, 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

Re-arranging atomic orbitals.

Figure 1: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}, χ3sF\chi_{3s}^{F} \chi_{1s}^{H}, χ2sH\chi_{2s}^{H} \chi_{1py}^{F}, χ2pyF\chi_{2py}^{F}, \chi_{1pz}^{F}, χ2pzF\chi_{2pz}^{F}, \chi_{1px}^{F}, χ2pxF\chi_{2px}^{F}

where orbitals classified as occupied are high-lighted in red. The T0T_0 transformation re-arranges the orbitals as follows:

\chi_{1s}^{F}

,

\chi_{2s}^{F}, χ3sF\chi_{3s}^{F}, \chi_{1px}^{F}

,

\chi_{1py}^{F}

,

\chi_{1pz}^{F}, χ2pxF\chi_{2px}^{F}, χ2pyF\chi_{2py}^{F},χ2pzF\chi_{2pz}^{F}, \chi_{1s}^{H}, χ2sH\chi_{2s}^{H}

rearranged_indices = vlx.aoindices.get_basis_function_indices_of_atoms(molecule, basis)

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

Diagonalization of the atomic blocks.

Figure 2: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

S0=LLTS_0 = L L^T

where LL is a lower triangular matrix. We get

I=L−1S0[LT]−1I = L^{-1} S_0 [L^T]^{-1}

and identify

T1=[LT]−1T_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

Diagonalization of the occupied and virtual blocks.

Figure 3: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

Projecting out the occupied space from the virtual.

Figure 4: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

Diagonalization the virtual subspace.

Figure 5: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

⟨Ω^⟩=tr[ΩD]\langle \hat{\Omega} \rangle = \mathrm{tr} \big[ \Omega D \big]
(3)

where Ω\Omega and DD 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

∣ϕ⟩=∣χ⟩C=∣χL⟩T−1C| \boldsymbol{\phi} \rangle = | \boldsymbol{\chi} \rangle C = | \boldsymbol{\chi}^\mathrm{L} \rangle T^{-1} C
(4)

we identify

CL=T−1CC^\mathrm{L} = T^{-1} C
(5)

and consequently

DL=T−1D [T†]−1=T−1D [T−1]†D^\mathrm{L} = T^{-1} D \, [T^\dagger]^{-1} = T^{-1} D \, [T^{-1}]^\dagger
(6)

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

ΩL=T†ΩT\Omega^\mathrm{L} = T^\dagger \Omega T
(7)

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

⟨Ω^⟩=∑AΩAL;ΩAL=∑BΩABL;ΩABL=∑α∈A∑β∈BΩαβLDβαL\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}
(8)

where AA and BB 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 tr[SD]Q = -e \, \mathrm{tr} \big[ S D \big]
(9)

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

QA=−e∑α∈ADααLQ_A = - e \sum_{\alpha \in A} D^\mathrm{L}_{\alpha\alpha}
(10)
# atomic numbers and therefore nuclear charge in atomic units
Z_A = molecule.get_element_ids()
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 Ω^=−e z^\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

⟨μ^⟩=tr[μD]\langle \hat{\mu} \rangle = \mathrm{tr} \big[ \mu D \big]
(11)
dipole_mats = vlx.compute_electric_dipole_integrals(molecule, basis)

mu = dipole_mats[2]  # mu_z

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.

Bond midpoint coordinates.

Figure 6:Bond midpoint coordinates.

We introduce

RAB=RA+RB2\boldsymbol{R}_{AB} = \frac{\boldsymbol{R}_{A} + \boldsymbol{R}_{B}}{2}
(12)

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

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

coords = molecule.get_coordinates_in_bohr()

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

μ(AB)L=−e∑α∈A∑β∈B⟨χαL∣(r^−RAB)∣χβL⟩DβαL\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}
(13)

It can be determined from the relation

μ(AB)L=μABL−QARAδAB\boldsymbol{\mu}^\mathrm{L}_{(AB)} = \boldsymbol{\mu}^\mathrm{L}_{AB} - Q_A \boldsymbol{R}_{A} \delta_{AB}
(14)
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

α(−ω;ω)=−⟨⟨μ^;μ^⟩⟩ω=−μ[1]†N(ω)\alpha(-\omega;\omega) = - \langle\langle \hat{\mu}; \hat{\mu} \rangle\rangle_\omega = -{\mu^{[1]}}^\dagger N(\omega)
(15)

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

The strucure of the property gradient is

μ[1]=(g−g∗)\mu^{[1]} = \begin{pmatrix} g \\ -g^* \\ \end{pmatrix}
(16)

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

E^si†=12(a^s,α†a^i,α+a^s,β†a^i,β)\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)
(17)

we get

gsi=⟨0∣[E^si,μ^]∣0⟩=2μsig_{si} = \langle 0 | \big[ \hat{E}_{si}, \hat{\mu} \big] | 0 \rangle = \sqrt{2} \mu_{si}
(18)

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

N(ω)=(E[2]−ℏωS[2])−1μ[1]N(\omega) = \Big( E^{[2]} - \hbar \omega S^{[2]} \Big)^{-1} \mu^{[1]}
(19)

The structure of the response vectors is

N(ω)=(Z−Y∗)N(\omega) = \begin{pmatrix} Z \\ -Y^* \\ \end{pmatrix}
(20)

and in the static limit, we have Y=ZY = 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

κ=(0ZY∗0)\kappa = \begin{pmatrix} 0 & Z \\ Y^* & 0 \\ \end{pmatrix}
(21)
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

α(−ω;ω)=tr[μ κ]=∑α,βμαβDαβ(1)\alpha(-\omega; \omega) = \mathrm{tr}\big[ \mu \,\kappa \big] = \sum_{\alpha,\beta} \mu_{\alpha\beta} D^{(1)}_{\alpha\beta}
(22)

where we have introduced the perturbed density matrix

Dαβ(1)=∑p,qCαp∗CβqκqpD^{(1)}_{\alpha\beta} = \sum_{p,q} C^*_{\alpha p} C_{\beta q} \kappa_{qp}
(23)

It can also be exressed in the LoProp basis

α(−ω;ω)=∑A,BαABL;αABL=∑α∈A∑β∈BμαβLDβαL,(1)\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}
(24)
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

α(AB)L=−e∑α∈A∑β∈B⟨χαL∣(r^−RAB)∣χβL⟩DβαL,(1)\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}
(25)

It can be determined from the relation

α(AB)L=αABL−ΔQARAδAB\alpha^\mathrm{L}_{(AB)} = \alpha^\mathrm{L}_{AB} - \Delta Q_A \boldsymbol{R}_{A} \delta_{AB}
(26)

where the external field derivative of the change in charge on atom AA is

ΔQA=−e∑α∈ADααL,(1)\Delta Q_A = -e \sum_{\alpha \in A} D^{\mathrm{L},(1)}_{\alpha\alpha}
(27)

As charge is preserved, we have

∑AΔQA=0\sum_A \Delta Q_A = 0
(28)
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, ΔQAB\Delta Q_{AB}, as the external field derivative of the charge transfer from center AA to BB. Although this quantity is not uniquely defined, it has to fulfill

∑A≠BΔQAB=ΔQA\sum_{A \neq B} \Delta Q_{AB} = \Delta Q_{A}
(29)

and

ΔQBA=−ΔQAB\Delta Q_{BA} = - \Delta Q_{AB}
(30)

In the LoProp scheme, it assumed that ΔQAB\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

L(ΔQAB,λA)=∑AB[ΔQAB]2f(ΔRAB)+∑AλA(∑A≠BΔQAB−ΔQA)\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)
(31)

where ΔRAB=∣RA−RB∣\Delta R_{AB} = |\boldsymbol{R}_A - \boldsymbol{R}_B| and ff is a penatly function preventing long-distance charge transfer chosen as

f(r)=exp⁡(α(rRABS+RBBS)2)f(r) = \exp\Big(\alpha \Big(\frac{r}{R_A^\mathrm{BS} + R_B^\mathrm{BS}}\Big)^2\Big)
(32)

The constant α\alpha is chosen equal to 2.0 and RABSR_A^\mathrm{BS} is the Bragg–Slater radius of atom AA.

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λ=ΔQL \lambda = \Delta Q
(33)

where

LAB=12f(ΔRAB)+CL_{AB} = \frac{1}{2 f(\Delta R_{AB})} + C
(34)

and

LAA=−∑A≠B12f(ΔRAB)+CL_{AA} = -\sum_{A \neq B} \frac{1}{2 f(\Delta R_{AB})} + C
(35)

The constant CC is determined from the requirement of charge conservation, or zero net charge transfer.

Once the Lagrange multipliers are determined, we have

ΔQAB=−λA−λB2f(ΔRAB)\Delta Q_{AB} = - \frac{\lambda_{A} - \lambda_{B}}{2 f(\Delta R_{AB})}
(36)
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 zzzz-component of the bond polarizability is defined by

αABbond=ΔQABΔRAB,z\alpha^\mathrm{bond}_{AB} = \Delta Q_{AB} \Delta \boldsymbol{R}_{AB, z}
(37)
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 AA is obtained as

αAL=∑B(α(AB)L+12αABbond)\alpha^\mathrm{L}_A = \sum_B \Big( \alpha^\mathrm{L}_{(AB)} + \frac{1}{2} \alpha^\mathrm{bond}_{AB} \Big)
(38)

and it fulfills

∑AαAL=α(−ω;ω)\sum_A \alpha^\mathrm{L}_A = \alpha(-\omega;\omega)
(39)
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.
References
  1. Gagliardi, L., Lindh, R., & Karlström, G. (2004). Local properties of quantum chemical systems: The LoProp approach. J. Chem. Phys., 121, 4494–4500. 10.1063/1.1778131
Environment
RESP charges
Environment
Embedding potentials