<jupyter_start><jupyter_text>Exploring molecular interactions using electronic structure methods. For today's activity we will use the package of programs psi4, so we will need to install it first. The psi4 python module can be installed using the anaconda python platform. To do this download anaconda in:https://www.anaconda.com/download/linuxYou can download either the 2.7 or 3.6 verison depending on your preference. Once downloaded and installed it is recomended to create a psi4 environment form whichthe psi4 module can be loaded. Run the comand:**conda create -n p4env python=2.7 psi4 psi4-rt -c psi4/label/dev -c psi4**This will install all the psi4 binaries and a python module which can be imported from the notebook.To finish the installation you need to provide a scratch directory in your .bashrc, for example:**export PSI_SCRATCH=/home/user_name/scratch/psi4**Finally source your .bashrc and activate the environment**source activate p4env**Now you may open jupyter-notebook and install psi4.<jupyter_code>import psi4 
import numpy as np<jupyter_output><empty_output><jupyter_text>1. Compute the energy of a diatomic molecules.As a first example we will compute the scf energy of the diatomic molecule hydrogen fluoride (HF):<jupyter_code># ==> Opciones Básicas Psi4 <==
# Memoria
psi4.set_memory(int(5e8))
numpy_memory = 500

# Output
psi4.core.set_output_file('output.dat', False)

# Geometría
hf_mol = psi4.geometry("""
0 1
H
F 1 0.917
""")

#psi4.optimize('scf/cc-pvtz')
#di_mat = hf_mol.distance_matrix()
#print np.asarray(di_mat)
energy_hf_mol , wfn_hf_mol = psi4.energy('mp2/cc-pvtz', return_wfn=True)
print energy_hf_mol<jupyter_output><empty_output><jupyter_text>This corresponds to the Born-Oppenheimer scf energy of HF. Now it's your turn. Compute the single point energy of F$_2$ and N$_2$ using the cc-pVDZ and cc-pVTZ.   2. Compute the dipole and quadrupole moment of diatomic molecules Since we are interested in studying long range  molecular interactions using classical electrodynamics, it is necessary to compute the dipole and quadrupole moments. Quantum mechanically the dipole can be computed using the one electron dipole operator:\begin{equation}\hat{\mu} = \sum_i q_i r_i \end{equation}where $q_i$ is the charge of the particle and $r_i$ is the position vector of the particle.  The dipole moment can be computed using the wavefunction through the expectation value of the operator $\mu$.\begin{equation}\mu =  \end{equation}In psi4 we can obtain the dipole moment from the wavefunction object that was defined above<jupyter_code>psi4.oeprop(wfn_hf_mol, 'DIPOLE', 'QUADRUPOLE', title='HF SCF')

mux = psi4.core.get_variable('HF SCF DIPOLE X') # in debye
muy = psi4.core.get_variable('HF SCF DIPOLE Y')
muz = psi4.core.get_variable('HF SCF DIPOLE Z')
quad_zz = psi4.core.get_variable('HF SCF QUADRUPOLE ZZ')
print muz
print quad_zz
mu = (np.sqrt(mux**2 + muy**2 + muz**2))
print mu<jupyter_output><empty_output><jupyter_text>3.  Compute a potential energy surface of HF dimer.In order to study the physical interactions between two molecules it is convenient to draw a potential energy surface along the interaction coordinate. In this section we will obtain a potential energy profile for the most favorable dipole-dipole interaction, which is the horizontal orientation with opposing dipole vectors, HF---HF. First we need to define a list containing the distances between both dimers for which the energy will be obtained.<jupyter_code>hf_dimer = psi4.geometry("""
        0 1
        H
        F 1 0.917
        H 2 R 1 180.0
        F 3 0.917 2 180.0 1 0.0
        """)<jupyter_output><empty_output><jupyter_text>Next, we write a loop and in each step of the loop we compute the energy at the mp4 level of theory.<jupyter_code>energy = []
dist   = []

Rval = np.arange(1.5,10.0,0.1)

for d in Rval:
    hf_dimer.R = d
    psi4.set_options({'freeze_core': 'True'})
    en = psi4.energy('scf/cc-pvtz')
    print en
    print d
    energy.append(en)
    dist.append(d)<jupyter_output><empty_output><jupyter_text>Now we are ready to plot the potential energy profile. We will use the matplotlib python library for this purpose. The function ref_cero_kcal transforms the energy which is in hartee to kcal/mol and takes the energy of the dimer with the farthest separation as the reference energy.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline
def ref_cero_kcal(en_list):
    energy_kcal = []
    for x in range(len(en_list)):
        energy_kcal.append((en_list[x] - en_list[-1])*627.51)
    return energy_kcal


energy_kcal = ref_cero_kcal(energy)    
plt.plot(dist,energy_kcal)<jupyter_output><empty_output><jupyter_text>**Exercise:** Using the data from the plot, make an estimate of the interaction energy in kcal/mol. Corroborate your results computing the interaction energy using the super-molecule approach. ( Hint: 1 hartree = 627.51 kcal/mol). What percentage of the total energy is the intermolecular interaction energy?<jupyter_code>hf_dimer_min = psi4.geometry("""
        0 1
        H
        F 1 0.917
        H 2 1.980 1 180.0
        F 3 0.917 2 180.0 1 0.0
        """)

energy_hf_dimer_min , wfn_hf_dimer_min = psi4.energy('mp2/cc-pvtz', return_wfn=True)
print energy_hf_dimer_min

int_energy = (energy_hf_dimer_min - 2*energy_hf_mol)

print int_energy*627.51

print("Percentage: "+str((int_energy/energy_hf_dimer_min)*100))<jupyter_output><empty_output><jupyter_text>4. The effect of electron correlationSince we know that the dispersion energy arises only when the motion of the electrons on the different molecular fragments are correlated, it is possible to evalute the importance of disperison efects by usingmethods that include diferent levels of electron correlation. **Exercices**: Construct a potential energy curve at the SCF, MP2 and MP4 level of theory and asses the importance of including electron correlation along the intermolecular coordinate.<jupyter_code>energy_scf = []
energy_mp2 = []
energy_mp4 = []
dist = []

for d in Rval:
    hf_dimer.R = d
    psi4.set_options({'freeze_core': 'True'})
    en_scf = psi4.energy('scf/cc-pvtz')
    en_mp2 = psi4.energy('mp2/cc-pvtz')
    en_mp4 = psi4.energy('mp4/cc-pvtz')
    print en_scf
    print d
    energy_scf.append(en_scf)
    energy_mp2.append(en_mp2)
    energy_mp4.append(en_mp4)
    dist.append(d)
energy_scf_kcal = ref_cero_kcal(energy_scf)  
energy_mp2_kcal = ref_cero_kcal(energy_mp2)  
energy_mp4_kcal = ref_cero_kcal(energy_mp4)  

plt.plot(dist,energy_scf_kcal)
plt.plot(dist,energy_mp2_kcal)
plt.plot(dist,energy_mp4_kcal)<jupyter_output><empty_output><jupyter_text>5. Energy decomposition with SAPTSo far we have only computed the interaction energy as a whole, however we know that we candecompose the interaction energy into contributions stemming from different physical phenomena, in particular, electrostatic, induction, dispersion and exchange interactions. First we will run a SAPT0 calculation to obtain the total interaction energy and compare it to the results we obtained above.<jupyter_code>hf_dimer_sapt = psi4.geometry('''
     0 1
     H  0.0 0.0 0.0
     F  0.0 0.0 0.917
     --
    0 1
     H  0.0 0.0 2.951
     F  0.0 0.0 3.868
     symmetry c1
''')

psi4.set_options({'basis': 'jun-cc-pVDZ',
                  'e_convergence': 1e-8,
                  'd_convergence': 1e-8})

psi4.energy('sapt0')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08289455920075604, 'fraction_numerical': 0.028621574186580262, 'mean_word_length': 4.427106227106227, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 38, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '483298', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2591, 'n_tokens_neox': 2421, 'n_words': 828, 'path': 'Theory_electronic_structure_day3/m_inter.ipynb', 'repo_name': 'QCMM/workshop2017', 'size': '7407'}