Script to get mean atomic contributions to heat-of-formation

lsq.py

#!/usr/bin/env python2
import ezpickle
import numpy as np

Q  = dict()
Q ["H"]  = 0
Q ["C"]  = 1
Q ["N"]  = 2
Q ["O"]  = 3
Q ["S"]  = 4
#Q ["Si"] = 5
# Q ["Ge"] = 6

def get_mean_atomic_contribution(mols, feature):

    feature = np.array(feature)

    nm = len(mols)
    nq = len(Q)

    # Make a matrix of number of atoms for each molecule
    at = np.zeros((nm,nq))

    for i, mol in enumerate(mols):
        # print i, mol.atomtypes
        for j, atomtype in enumerate(mol.atomtypes):

            at[i, Q[atomtype]] += 1.0

    # Solve n_atoms x mean_energy_per_atom = energy_for_molecule
    feature_atomic = np.linalg.lstsq(at,feature)[0]

    # Expand into atomic-averaged feature
    feature_mean = np.dot(at, feature_atomic)

    feature_compensated = feature - feature_mean

    return feature_compensated, feature_atomic
    

if __name__ == "__main__":

    mols = ezpickle.load("../atomic_krr/mols.cpickle")    

    es = [mol.energy for mol in mols]

    print get_mean_atomic_contribution(mols, es)