Example for numerical alchemical derivatives with the FCHL representation in QML.

alchemy_example.py

#!/usr/bin/env python

from __future__ import print_function

import numpy as np

import qml
import qml.fchl

from qml.fchl import generate_representation
from qml.fchl import get_atomic_kernels






def get_overlap(Z, width):

    n = len(Z)
    S = np.zeros((n,n))

    for i, Zi in enumerate(Z):
        for j, Zj in enumerate(Z):

            S[i,j] = np.exp(-(Zi - Zj)**2 / (2*width**2))

    return S


    


if __name__ == "__main__":

    # Nuclear charge displacement
    dZ = 0.02

    # Nuclear charges and displacements for 10 first atoms in periodic system
    Z = sum([[i-dZ, float(i), i+dZ] for i in range(1, 10)], [])
    print(Z)

    np.set_printoptions(linewidth=666, suppress=True)

    # Large Gaussian overlap between adjacent elements
    width = 1.5
    S = get_overlap(Z, width)

    print(S[:12,:12])
    
    # Small Gaussian overlap between adjacent elements
    width = 1.5
    S = get_overlap(Z, width)

    print(S[:12,:12])

    # Water molecule
    # 3
    #
    # O 1.464 0.707 1.056
    # H 0.878 1.218 0.498
    # H 2.319 1.126 0.952

    coordinates = np.array([
        [1.464, 0.707, 1.056],
        [0.878, 1.218, 0.498],
        [2.319, 1.126, 0.952]
    ])

    # Index, starts from 1, not zero.
    nuclear_indices = [23, 2, 2]


    # Generate representation for the water molecule
    rep = generate_representation(coordinates, nuclear_indices)
    rep = rep[:3]

    # Kernel elements for the atoms in the water molecule, 
    # "[2.5]" is the kernel width for the gaussian kernel.
    # alchemy=S means that the alchemical overlap is read from the matrix, S
    K = get_atomic_kernels(rep, rep, [2.5], alchemy=S)
    print(K)

    # The indices of the nuclear charges displaced by +/- dZ
    nuclear_indices_minus = [24, 2, 2]
    nuclear_indices_plus  = [22, 2, 2]
   
    # The displaced representations
    rep_minus = generate_representation(coordinates, nuclear_indices_minus)[:3]
    rep_plus  = generate_representation(coordinates, nuclear_indices_plus)[:3]


    # Kernel elements with the displaced nuclear charges
    K_minus = get_atomic_kernels(rep_minus, rep, [2.5], alchemy=S)
    K_plus  = get_atomic_kernels(rep_plus, rep, [2.5], alchemy=S)

    print(K_minus)
    print(K_plus)

    # The resulting numerical derivative wrt. the nuclear charge of Oxygen
    dKdZ = (K_plus - K_minus) / (2 * dZ)

    print(dKdZ)