<jupyter_start><jupyter_text>Quantum Chemistry in Python===============In this lesson we will explore how to perform Quantum Chemistry computations from Python,as well as several data science applications within the field of quantum chemistry.  To do this, we will utilize the quantum chemistry package Psi4, which has a Python front end and API, as well as several Python libraries for data analysis, visualization, and machine learning.  Quantum Chemistry BasicsThe goal of the field of quantum chemistry is to model, predict, and describe chemical behaviorfrom first principles using quantum mechanics.  To do so exactly is computaionally impossiblefor all but the smallest atomic and molecular systems, due to the exponential cost of solvingthe exact electronic Schrodinger equation.  Therefore, in quantum chemistry, we employ approximatemethods which* are _computationally tractible_ for "real" chemical systems of interest and * provide an _accurate description of the chemical behavior_ to be studied.<jupyter_code># ==> Import libraries <==
import psi4                               # Quantum chemistry package 
import time                               # For timing code
import numpy as np                        # Numerical operations
import pandas as pd                       # Data handling
import seaborn as sns                     # Plotting & visualization
import matplotlib.pyplot as plt           # Plotting & visualization
import utilities.plotting as fancyplt     # Custom library for advanced plots
import utilities.misc as misc             # Custom functions for stats, IEs, etc.

idx = pd.IndexSlice                       # Slice pandas.DataFrame objects across multiindex

# Matplotlib plots within notebook
%matplotlib inline
# ==> Example 1: Construct He...He Potential Curve <==

# Define atomic coordinates for He dimer
He_dimer = """
0 1
He 0 0 0    # First He at the origin
--
He 0 0 {}   # Other He displaced along z-axis

units angstrom
"""

psivariables = {}

# Iterate over He...He separation distances in z axis
for r in np.arange(2.6, 10, 0.1):
    mol = psi4.geometry(He_dimer.format(r))                      # Create Psi4 molecule with displaced He dimer
    ie = psi4.energy('b3lyp-d3mbj/aug-cc-pvdz', bsse_type='cp')  # Compute CP-corrected interaction energy
    psivariables[r] = psi4.core.variables()                      # Save Psi4 internal variables

b3lypdf = pd.DataFrame(data=psivariables).T # Store Psi4 internal variables inside Pandas DataFrame
b3lypdf *= psi4.constants.hartree2kcalmol   # Convert energy units from atomic units (Hartree) to kcal/mol

# Plot potential curve directly from Pandas dataframe
plot_options = {'style': 'bo-',
                'title': 'He$_2$ Interaction Energies',
                'legend': True
               }

b3lypdf['CP-CORRECTED 2-BODY INTERACTION ENERGY'].plot(**plot_options)
plt.xlabel('Interatomic Separation ($\AA$)')
plt.ylabel('Interaction Energy (kcal/mol)')

# Save dataframe to disk for future use
b3lypdf.to_pickle('data/He2-b3lyp-d3mbj-curve.pkl')
# ==> Plot IE curve directly from Pandas <==
plot_options = {'style': 'bo-',
                'title': 'He$_2$ Interaction Energies',
                'legend': True
               }

b3lypdf['CP-CORRECTED 2-BODY INTERACTION ENERGY'].plot(**plot_options)
plt.xlabel('Interatomic Separation ($\AA$)')
plt.ylabel('Interaction Energy (kcal/mol)')<jupyter_output><empty_output><jupyter_text>Improving QC Model Accuracy using Transfer Learning Symmetry-Adapted Perturbation Theory: Physically Motivated Interaction Energies from First PrinciplesSymmetry-Adapted Perturbation Theory (SAPT) can provide a more chemically intuitive description of intermolecular interaction energies by providing a breakdown of the IE into phisically meaningful components:$$ IE_{\rm SAPT0} = E_{\rm elst} + E_{\rm exch} + E_{\rm ind} + E_{\rm disp} $$To illustrate this method, we can visualize the total IE and components of our Helium dimer using SAPT:<jupyter_code># ==> Example 02: SAPT IE & Components for He2 <==

# Compute SAPT IE for He2 along potential curve
for r in np.arange(2.6, 10, 0.1):
    mol = psi4.geometry(He_dimer.format(r))                      # Create Psi4 molecule with displaced He dimer
    ie = psi4.energy('sapt0/jun-cc-pvdz', molecule=mol)          # Compute SAPT IE
    psivariables[r] = psi4.core.variables()                      # Save Psi4 internal variables

saptdf = pd.DataFrame(data=psivariables).T      # Store Psi4 internal variables inside Pandas DataFrame
saptdf *= psi4.constants.hartree2kcalmol        # Convert energy units from atomic units (Hartree) to kcal/mol

# Visualize SAPT components along potential scan
terms = {'SAPT ELST ENERGY': 'ro-',
         'SAPT EXCH ENERGY': 'go-',
         'SAPT IND ENERGY': 'bo-',
         'SAPT DISP ENERGY': 'mo-',
         'SAPT TOTAL ENERGY': 'ko-'
        }

saptdf[terms.keys()].plot(style=list(terms.values()), markersize=2)
plt.xlabel('Interatomic Separation ($\AA$)')
plt.ylabel('Interaction Energy (kcal/mol)')
plt.title('SAPT0 IE & Components')

# Plot SAPT IE vs. B3LYP-D3M(BJ) IE
df = pd.concat([saptdf['SAPT TOTAL ENERGY'], b3lypdf['CP-CORRECTED 2-BODY INTERACTION ENERGY']], axis=1)
df.columns = ['SAPT0', 'B3LYP-D3M(BJ)']
df.plot(style='o-')
plt.xlabel('Interatomic Separation ($\AA$)')
plt.ylabel('Interaction Energy (kcal/mol)')
plt.title('SAPT0 vs. B3LYP-D3M(BJ) IE')

# Save SAPT data for later 
saptdf.to_pickle('data/He2-sapt0-curve.pkl')
# ==> Visualize the components of SAPT along potential scan <==

terms = {'SAPT ELST ENERGY': 'ro-',
         'SAPT EXCH ENERGY': 'go-',
         'SAPT IND ENERGY': 'bo-',
         'SAPT DISP ENERGY': 'mo-',
         'SAPT TOTAL ENERGY': 'ko-'
        }

saptdf[terms.keys()].plot(style=list(terms.values()), markersize=2)
plt.xlabel('Interatomic Separation ($\AA$)')
plt.ylabel('Interaction Energy (kcal/mol)')
plt.title('SAPT0 IE & Components')

# Plot SAPT IE vs. B3LYP-D3M(BJ) IE
df = pd.concat([saptdf['SAPT TOTAL ENERGY'], b3lypdf['CP-CORRECTED 2-BODY INTERACTION ENERGY']], axis=1)
df.columns = ['SAPT0', 'B3LYP-D3M(BJ)']
df.plot(style='o-')
plt.xlabel('Interatomic Separation ($\AA$)')
plt.ylabel('Interaction Energy (kcal/mol)')
plt.title('SAPT0 vs. B3LYP-D3M(BJ) IE')<jupyter_output><empty_output><jupyter_text>SAPT0-D: A More Efficient SAPT ApproachFor large molecules the SAPT approach can become prohibitively expensive, because the cost of computing the SAPT energy scales as $N^5$, where $N$ is proportional to the total number of atoms. To illustrate this, consider the following set of computations computing the interaction between two growing chains of He atoms:<jupyter_code># ==> Visualize N^5 SAPT Scaling <==
chainA = ['He 0 0 0\n', 'He 0 0 2\n', 'He 0 0 4\n', 'He 0 0 6\n', 'He 0 0 8\n', 'He 0 0 10\n']
chainB = ['He 0 3 0\n', 'He 0 3 3\n', 'He 0 3 4\n', 'He 0 3 6\n', 'He 0 3 8\n', 'He 0 3 10\n']

molstr = """
0 1
{}
--
0 1
{}

units angstrom
"""

monoA = ''
monoB = ''
times = []
natoms = []
for i in range(len(chainA)):
    monoA += chainA[i]
    monoB += chainB[i]
    t0 = time.time()
    psi4.energy('sapt0/jun-cc-pvdz', mol=psi4.geometry(molstr.format(monoA, monoB)))
    times.append(time.time() - t0)
    natoms.append(2 * i + 1)
    
# Plot time versus number of atoms, fit to N^5 curve
f = np.poly1d(np.polyfit(natoms, times, deg=5))
plt.plot(natoms, f(natoms), label='${\cal O}(N^5)$')
plt.plot(natoms, times, 'ro', label='SAPT0 Time', zorder=10)
plt.title('Scaling of SAPT0')
plt.xlabel('Number of Atoms')
plt.ylabel('Time (s)')
plt.legend()<jupyter_output><empty_output><jupyter_text>The term responsible for the ${\cal O}(N^5)$ scaling of SAPT0 is the dispersion term,$E_{\rm disp}$, whereas all other terms scale as ${\cal O}(N^4)$. We could therefore seemingly reduce the computational scaling of SAPT0 by an order of magnitude by simplyneglecting the dispersion term.  But how could such an approach still produce accurateinteraction energies?We may be in luck! There exist empirical "dispersion corrections" which have been proposedto approximately add long-range electron correlation (the physical phenomenon responsiblefor the dispersion interaction) to mean-field quantum chemical methods like Density FunctionalTheory which are (either partially or entirely) neglecting this effect.  Importantly, methodswhich employ these dispersion corrections satisfy both of our desired properties in quantumchemistry: in general, the dispersion-corrected method acheives higher accuracy in the same amount of time as the uncorrected version of the method. This increased accuracy is achievedby training the empirical model to reproduce high-quality references, computed using methodswhich are _significantly_ more expensive than the model being corrected.  Most often, thesereference energies are computed using the quantum chemical "gold standard," CCSD(T). So, instead of simply abandoning the expensive SAPT0 dispersion, let's replace it with an empirical dispersion correction that can be parameterized to reproduce very high accuracyinteraction energies:  $$ IE_{\rm SAPT0-D} = E_{\rm elst} + E_{\rm exch} + E_{\rm ind} + E^{\rm -D}_{\rm disp} $$We'll refer to this method as SAPT0-D, for "dispersion-corrected SAPT0".<jupyter_code># Import SAPT, SAPT-D, & CCSD(T) Data
saptd = pd.read_pickle('data/saptd-full.pkl')

print(saptd.shape)
saptd.head()
# ==> Plotting Signed Errors for SAPT0 & SAPT0-D vs. CCSD(T) <==

rawdata=pd.concat([saptd.loc[:,idx['D3M',['SAPT0 Error','Benchmark']]], saptd.loc[:,idx[:,['SAPT0-Zero3 Error', 'SAPT0-BJ Error', 'SAPT0-TT Error']]]], axis=1)

# Clean data for potentially missing (or too small) benchmarks
rawdata = rawdata[abs(rawdata.loc[:,idx['D3M','Benchmark']]) >= 1.e-6] 
data=rawdata

se = pd.DataFrame(index=data.index, columns=['SAPT0','SAPT0-D3M','SAPT0-D3BJ'])#,'SAPT0-TT'])
se['SAPT0'] = data['D3M','SAPT0 Error']
se['SAPT0-D3M'] = data['D3M','SAPT0-Zero3 Error']
se['SAPT0-D3BJ'] = data['D3BJ','SAPT0-BJ Error']

bounds = [.59, 0.59, 0.4, 0.4]
insetopts = {'keep': 95,
             'fliers': False}
options={
            'labels': {'y axis': 'Signed Error (kcal/mol)', 'title': f'Quality of SAPT0 & SAPT0-D IEs'},
        }
fancyplt.violin(se, options)
# ==> Ternaries for Signed Errors <==
dispmodels = {'D3M': 'SAPT0-Zero3 Error', 'D3BJ': 'SAPT0-BJ Error' , 'TT': 'SAPT0-TT Error'}
# SAPT0
data = pd.DataFrame(index=saptd.index, columns=['Elst', 'Ind', 'Disp', 'Error'])
data['Error'] = saptd['D3M', 'SAPT0 Error']
data['Elst'] = saptd['D3M', 'Elst']
data['Ind'] = saptd['D3M', 'Ind']
data['Disp'] = saptd['D3M', 'Disp']
filtered = misc.qfilter(data, filter_on=['Error'], q=[0.05,0.995], ret='outlier', keep_index=True)
fancyplt.ternary(filtered.dropna(subset=['Error']), colors='error', colorbar='symmetric', title='SAPT0 SE')

# SAPT0-D3(0)
data = pd.DataFrame(index=saptd.index, columns=['Elst', 'Ind', 'Disp', 'Error'])
data['Error'] = saptd['D3M', 'SAPT0-Zero3 Error']
data['Elst'] = saptd['D3M', 'Elst']
data['Ind'] = saptd['D3M', 'Ind']
data['Disp'] = saptd['D3M', 'Disp']
filtered = misc.qfilter(data, filter_on=['Error'], q=[0.05,0.995], ret='outlier', keep_index=True)
fancyplt.ternary(filtered.dropna(subset=['Error']), colors='error', colorbar='symmetric', title='SAPT0-D3M(0) SE')

# SAPT0-D3(BJ)
data = pd.DataFrame(index=saptd.index, columns=['Elst', 'Ind', 'Disp', 'Error'])
data['Error'] = saptd['D3BJ', 'SAPT0-BJ Error']
data['Elst'] = saptd['D3M', 'Elst']
data['Ind'] = saptd['D3M', 'Ind']
data['Disp'] = saptd['D3M', 'Disp']
filtered = misc.qfilter(data, filter_on=['Error'], q=[0.05,0.995], ret='outlier', keep_index=True)
fancyplt.ternary(filtered.dropna(subset=['Error']), colors='error', colorbar='symmetric', title='SAPT0-D3M(BJ) SE')

# SAPT0-TT
data = pd.DataFrame(index=saptd.index, columns=['Elst', 'Ind', 'Disp', 'Error'])
data['Error'] = saptd['TT', 'SAPT0-TT Error']
data['Elst'] = saptd['D3M', 'Elst']
data['Ind'] = saptd['D3M', 'Ind']
data['Disp'] = saptd['D3M', 'Disp']
filtered = misc.qfilter(data, filter_on=['Error'], q=[0.05,0.995], ret='outlier', keep_index=True)
fancyplt.ternary(filtered.dropna(subset=['Error']), colors='error', colorbar='symmetric', title='SAPT0-TT SE')<jupyter_output><empty_output><jupyter_text>SAPT0-D(ML): Transfer Learning SAPT-D towards Quantum Chemical "Gold Standard"SAPT0-D successfully reduced the computational cost of SAPT0 from ${\cal O}(N^5)\rightarrow{\cal O}(N^4)$,however its accuracy with respect to CCSD(T) was not appreciably better than SAPT0.  In an attempt toincrease the accuracy of SAPT0-D without significantly increasing its complexity, we can employ "transferlearning," a concept which seeks to apply additional knowledge to improve an existing model, instead oftraining a completely new model:![blargh](media/transfer-learning.png)We can apply transfer learning to increase the accuracy of SAPT0-D by learning the remaining differencebetween SAPT0-D and CCSD(T), and then predicting this correction using a simple regression model.<jupyter_code># ==> Fitting Delta CCSD(T) with KRR <==
import keras
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import KFold
from sklearn.kernel_ridge import KernelRidge

seed=314
np.random.seed(seed)

# Read in SAPT-D data, reorganize slightly
data = pd.read_pickle("data/full-disponly-dashD.pkl")
cols = pd.MultiIndex.from_product([['D3BJ', 'D3M', 'TT'], ['BJ', 'Benchmark', 'Disp', 'Elst', 'Exch', 'Fitset', 'Ind', 'Reference', 'SAPT0', 'SAPT0 Error', 'SAPT0 UE', 'SAPT0 URE', 'SAPT0 WURE', 'SAPT0-BJ', 'SAPT0-BJ Error', 'SAPT0-BJ UE', 'SAPT0-BJ URE', 'SAPT0-BJ WURE', 'SAPT0-TT', 'SAPT0-TT Error', 'SAPT0-TT UE', 'SAPT0-TT URE', 'SAPT0-TT WURE', 'SAPT0-Zero3', 'SAPT0-Zero3 Error', 'SAPT0-Zero3 UE', 'SAPT0-Zero3 URE', 'SAPT0-Zero3 WURE', 'SAPTHF', 'TT', 'Zero3', 'dashD(ML)-cc', 'dashD(ML)-delta']])
saptdml = pd.DataFrame(columns=cols, index=data.index)

# Iterate over damping functions, fit separate \Delta CC
for damp in ['D3BJ', 'TT', 'D3M']:
    if damp == 'D3BJ':
        saptname = 'BJ'
    elif damp == 'D3M':
        saptname = 'Zero3'
    else: saptname=damp
    
    # Select features to train KRR on
    feature_columns = [(damp, f'SAPT0-{saptname}'), (damp, 'SAPTHF'), (damp, 'Elst'), (damp, 'Ind'),(damp, 'Exch'), (damp, saptname)]
    
    # Choose fitting set for KRR, use same as for SAPT-D
    fitset = data.loc[:,(damp, 'Fitset')].values
    cc_tot = data.loc[:,[(damp,'Reference')]].values[:,0] #USE FOR BAD EMPIRICAL DISP CORRECTION
    
    sapt0_d3_tot = data.loc[:, [(damp,f'SAPT0-{saptname}')]].values[:,0]
    inv_disp_tot = data.loc[:, [(damp,'Elst')]].values[:,0] / data.loc[:, [(damp,saptname)]].values[:,0]
    
    x_data = data.loc[:,feature_columns].values
    y_data = cc_tot - np.array(data.loc[:,feature_columns[0]].values)
    print(f"SAPT0-{damp}(ML) INPUT PARAMETERS :\n   {feature_columns}\n")
    
    x_tot = x_data
    y_tot = y_data
    
    x = []
    x_holdout = []
    y = []
    y_holdout = []
    
    inv_disp = []
    inv_disp_holdout = []
    inv_disp_data = []
    
    cc = []
    cc_holdout = []
    cc_data = []
    sapt0_d3 = []
    sapt0_d3_holdout = []
    sapt0_d3_data = []
    
    
    for i in range(len(y_tot)):
        inv_disp_data.append(inv_disp_tot[i])
        cc_data.append(cc_tot[i])
        sapt0_d3_data.append(sapt0_d3_tot[i])
        if fitset[i] == True:
            x.append(x_tot[i])
            y.append(y_tot[i])
            cc.append(cc_tot[i])
            sapt0_d3.append(sapt0_d3_tot[i])
            inv_disp.append(inv_disp_tot[i])
        elif fitset[i] == False:
            x_holdout.append(x_tot[i])
            y_holdout.append(y_tot[i])
            cc_holdout.append(cc_tot[i])
            sapt0_d3_holdout.append(sapt0_d3_tot[i])
            inv_disp_holdout.append(inv_disp_tot[i])
    
    # Perform K-fold cross validation for KRR hyperparameter fitting
    x = np.array(x)
    y = np.array(y)
    x_holdout = np.array(x_holdout)
    y_holdout = np.array(y_holdout)
    inv_disp = np.array([inv_disp])
    inv_disp_holdout = np.array([inv_disp_holdout])
    inv_disp_data = np.array([inv_disp_data])
    cc = np.array(cc)
    cc_holdout = np.array(cc_holdout)
    cc_data = np.array(cc_data)
    sapt0_d3 = np.array(sapt0_d3)
    sapt0_d3_holdout = np.array(sapt0_d3_holdout)
    sapt0_d3_data = np.array(sapt0_d3_data)
    x = np.hstack((x, inv_disp.T))
    x_holdout = np.hstack((x_holdout, inv_disp_holdout.T))
    x_data = np.hstack((x_data, inv_disp_data.T))
 
    k_fold = KFold(n_splits=10, shuffle=True, random_state=seed)
    k_fold.get_n_splits(x)
    
    avg_y = np.mean(y)
    cvscores = []
    
    # Fit KRR model for delta CC
    model = KernelRidge(kernel='polynomial', degree=3)
    model.fit(x,y)
    
    # Predict delta CC
    preds = model.predict(x_data)
    
    # Add delta CC to SAPT0-D to predict CC
    cc_pred = preds + sapt0_d3_data
    
    model_abs_errs = np.abs(cc_pred - cc_data)
    model_errs = cc_pred - cc_data
    sapt_abs_errs = np.abs(sapt0_d3_data - cc_data)
    sapt_errs = sapt0_d3_data - cc_data
    
    print(f"   SAPT0/{damp} MAE on predicting CCSD(T) : {np.mean(sapt_abs_errs)}")
    print(f"   SAPT0/{damp} mean signed error on predicting CCSD(T) : {np.mean(sapt_errs)}")
    print(f"   KRR + SAPT0/{damp} MAE on predicting CCSD(T) : {np.mean(model_abs_errs)}")
    print(f"   KRR + SAPT0/{damp} mean signed error on predicting CCSD(T) : {np.mean(model_errs)}\n")
   
    # Save predicted values to dataframe
    saptdml[(damp,'dashD(ML)-delta')] = pd.Series(preds).values
    saptdml[(damp,'dashD(ML)-cc')] = pd.Series(cc_pred).values
    
    # Save errors in predicted values to dataframe
    saptdml.loc[:, idx[damp, f'SAPT0-{damp}(ML) SE']] = saptdml.loc[:,idx[damp,'dashD(ML)-cc']] - saptd.loc[:, idx[damp,'Benchmark']]
    
# Write dataframe to disk to save for later
saptdml.to_pickle("data/dashD_ML.pkl")
# ==> Plotting Signed Errors for SAPT0 & SAPT0-D vs. CCSD(T) <==

rawdata=pd.concat([saptd.loc[:,idx['D3M',['SAPT0 Error','Benchmark']]], saptdml.loc[:,idx[:,['SAPT0-D3M(ML) SE', 'SAPT0-D3BJ(ML) SE', 'SAPT0-TT(ML) SE']]]], axis=1)

# Clean data for potentially missing (or too small) benchmarks
rawdata = rawdata[abs(rawdata.loc[:,idx['D3M','Benchmark']]) >= 1.e-6] 
data=rawdata

se = pd.DataFrame(index=data.index, columns=['SAPT0','SAPT0-D3M','SAPT0-D3BJ'])#,'SAPT0-TT'])
se['SAPT0'] = data['D3M','SAPT0 Error']
se['SAPT0-D3M'] = data['D3M','SAPT0-D3M(ML) SE']
se['SAPT0-D3BJ'] = data['D3BJ','SAPT0-D3BJ(ML) SE']

bounds = [.59, 0.59, 0.4, 0.4]
insetopts = {'keep': 95,
             'fliers': False}
options={
            'labels': {'y axis': 'Signed Error (kcal/mol)', 'title': f'Quality of SAPT0 & SAPT0-D(ML) IEs'},
        }
fancyplt.violin(se, options)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14814814814814814, 'fraction_numerical': 0.021817794806300552, 'mean_word_length': 4.088816680205794, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '450807', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7188, 'n_tokens_neox': 6824, 'n_words': 1762, 'path': 'Day05/Sherrill_ML-for-Drug-Binding/quantum-chemistry-in-python.ipynb', 'repo_name': 'kxu08/Bootcamp2019', 'size': '18792'}