Exercises

Introduction

Some general remarks at the beginning: There will be an accompanying GitHub Repository to this course. In this Repository, you can find all provided input files and geometries necessary to work on the tasks. In addition, you can also find the manuals of the program packages TURBOMOLE7.6, COSMOtherm, VASP5.4, and CRYTSAL14. The PSI4 documentation is available online.

Most of the other scripts give some general options if started via

[program] -h

Important

Use only your /tmp1/$USER/ folder for all calculations!

Noncovalent Interactions and Organic Solids

Partitioning noncovalent interactions

Exercise 1.1

Partition the interactions of three weakly bound dimers (namely the water dimer, the argon dimer, and the uracil dimer) and identify the dominant binding motifs.

Approach

We will use the symmetry-adapted perturbation theory (SAPT) implemented in PSI4. For your ease, we will provide the equilibrium geometries in XMol and TM formats. Choose the appropriate one. The SAPT ansatz gives the first and second-order complexation energies based on either an HF or a DFT monomer approximation:

\[E_{SAPT} = E_{pol} + E_{exch} + E_{ind,resp} + E_{exch-ind,resp} + E_{disp} + E_{exch-disp}\]

1. Calculate each system’s HF-SAPT/aug-cc-pVQZ and AC-PBE0-SAPT/aug-cc-pVQZ (second order) binding energy at the equilibrium distance. For the DFT calculation, you need to asymptotically correct the PBE0 functional shifting the ionization potential to a reference value. For the uracil dimer, you should use a smaller aug-cc-pVTZ basis. Discuss the problems that can arise by using a smaller basis. Classify the different binding situations by separating the first and second-order energy contributions. Explain why the differences between the HF and DFT-based SAPT energies increase in the order argon, water, and uracil.

   Technical procedure

   For the DFT-SAPT calculation, you need the difference between the experimental and calculated ionization energies of the monomers. You can find experimental ionization energies in the NIST database. Modify the provided input files to match your system. Modify the HF input file to choose the appropriate SAPT method for the above equation. Start the program by invoking

   psi4 [options]
   

   Hint

   PSI4 uses input files to perform calculations. The default input file name is input.dat. You can change this via command line options. PSI4 also uses an output file output.dat by default. Additional outputs may still be found in the standard output of your system. If PSI4 returns an error message about insufficient memory, try to provide the sufficient amount of memory to the calculation with the keyword --memory xxGB. You can find further information about PSI4 and SAPT calculations in the PSI 4 documentation.

2. Calculate the HF-SAPT2/aug-cc-pVQZ potential energy surface for the argon dimer. Discuss the different first and second order contributions at the different distances and plot the total electrostatic, exchange, induction, and dispersion contributions as well as the total SAPT2 interaction energy with respect to the Ar⋅ ⋅ ⋅Ar distance. What characteristic distance dependence do you see for the first order exchange \(E^{(1)}_{exch}\) and the second order dispersion \(E^{(2)}_{disp}\)?

   Technical procedure

   The distance scan can easily be performed via an external bash script. To get their functional form, plot the distance dependence and fit an appropriate function to the first order exchange and second order dispersion contribution.

Supermolecular approaches

Exercise 1.2

Calculate reference binding energies for the water and the argon dimer using the provided relaxed geometries.

Approach

The reference energies are calculated in the supermolecular approach with TURBOMOLE. The binding energy of a two-fragment system \(E_{bind}\) is given by the energy differences between the single fragments and the complete system.

\[E_{bind} = E^{AB} - E^{A} - E^{B}\]

The extrapolation to the basis set limit can be done by assuming a certain functional form.

\[E_{SCF}^{(X)} = E_{SCF}^{(\infty)} + A \cdot e^{-\alpha \sqrt{X}}\]

\[E_{corr}^{(X)} = E_{corr}^{(\infty)} + A \cdot X^{-\beta}\]

The following exponents have been optimized according to the different basis sets.

	\(\alpha_{2,3}\)	\(\beta_{2,3}\)	\(\alpha_{3,4}\)	\(\beta_{3,4}\)
cc-pVXZ	4.42	2.46	5.46	3.05
aug-cc-pVXZ	4.30	2.51	5.79	3.05
pc-X	7.02	2.01	9.78	4.09
def2-XZVP	10.39	2.40	7.88	2.97

1. Calculate the CCSD(T) energy (coupled cluster including singles, doubles and perturbative triples) in the aug-cc-pVTZ basis set. Why is the usage of augmented dunning basis sets reasonable?

2. Calculate the MP2 energy (Møller-Plesset second order perturbation theory) in the aug-cc-pVTZ and aug-cc-pVQZ basis sets.

   Technical procedure

   Use the TURBOMOLE input generator cefine_current to prepare the calculations (try -h to figure out which options you have to use). Generally add -sym c1 to your cefine call to avoid smmetry-related problems. First do a canonical Hartree-Fock calculation (ridft) and then compute the correlation energy (ccsdf12 for CCSD(T) and ricc2 for MP2).

3. Use the MP2 energies to extrapolate the CCSD(T) values to the complete basis set limit and discuss the results. Explain why this two-point extrapolation is justified. Compare the CCSD(T)/CBS(est.) energies with the plain Hartree-Fock values and the AC-PBE0-SAPT values. Which contributions to the binding energy are still missing, i.e. which errors are made?

Molecules in solution

Exercise 1.3

Calculate the equilibrium association free energy \(\Delta G_a\) of a molecular “tweezer” with tetracyanoquinone (TCNQ) at room temperature solvated in toluene.

Approach

The host-guest system is again treated in a supermolecular approach. The association free energy \(\Delta G_a\) is given by

\[\Delta G_{a} = \Delta E + \Delta G_{RRHO}^{T} + \Delta \delta G_{solv}^{T}(X)\]

with the electronic gas phase association energy \(\Delta E\), a correction to free energies in the rigid rotor harmonic oscillator approximation \(\Delta G_{RRHO}^{T}\), and a correction to the solvation free energy \(\Delta \delta G_{solv}^{T}(X)\). These contributions depend explicitly on the temperature and solvent.

Hint

Experimental value: \(\Delta G_{a}^{(298 K)} = -4.50\) kcal⋅mol^-1

1. Calculate the electronic energy contribution with the (two- and three-body dispersion) corrected meta-GGA density functional TPSS-D3^ATM(BJ) in the def-TZVP single particle basis set. Correct for the basis set superposition error via the geometrical counterpoise correction (gCP).

   Technical procedure

   TPSS-D3/def-TZVP geometries are provided. Calculate a TPSS single-point energy with TURBOMOLE. Calculate the D3^ATM(BJ) contribution with the dftd3 standalone program. (Attention: cefine sets a dispersion correction by default. Make sure that you don’t double-count it.) The counterpoise correction can be calculated with the gcp program. Use the -h (help) option of dftd3 and gcp to figure out the desired options. Make sure you are using the correct gcp version (v2.01), if not try loading the program with module load gcp/2.01

2. Calculate the vibrational contributions in the rigid rotor harmonic oscillator model at the semiempirical GFN2-xTB level.

   Technical procedure

   First, re-optimize the TPSS-D3 geometries at the GFN2-xTB level. Then calculate the second derivatives and read the thermodynamic functions printout in the RRHO approximation of the GFN2-xTB program. Further information is given in Section GFN-xTB.

3. Compute the solvent corrections with COSMO-RS.

   Technical procedure

   The COSMO-RS model can be used via TURBOMOLE and the cosmosolv_prak script (option -scf). The input file ~/.cosmothermrc has to be modified to match the correct solvent (see Section COSMOtherm).

Organic solids

Exercise 1.4

Calculate the sublimation enthalpy of the urea crystal at room temperature.

Approach

The sublimation enthalpy \(\Delta H_{sub}\) at temperature \(T\) is given by the enthalpy difference of the two phases:

\[\Delta H_{sub} = H^{g} - H^{s}\]

\[H^{i} = E_{el}^{i} + E_{trans}^{i} + E_{rot}^{i} + E_{vib}^{i} + pV\]

The gas phase is thereby denoted with the index \(g\), the solid phase with \(s\). The enthalpy \(H\) is the total internal energy \(E_{tot}\) with the inclusion of the volume work \(pV\). The indices correspond to the electronic (\(el\)), translational (\(trans\)), rotational (\(rot\)), and vibrational (\(vib\)) contributions to the total internal energy. The rotational and vibrational contributions to the gas phase free energy are typically modeled by an ideal gas. The vibrational contributions are treated in the harmonic approximation:

\[E_{vib}(T) = \sum_{k}^{modes} \left( \frac{\hbar\omega_k}{2} + \frac{\hbar\omega_k}{\text{exp}\left(\frac{\hbar\omega_k}{k_BT}\right)-1} \right)\]

Hint

Experimental value: \(\Delta H_{sub}^{(298 K)} = 22.42\) kcal⋅mol^-1

1. Convert the Crystallographic Information File (urea_113.cif) into CRYSTAL 14 fort.34 format.

   Technical procedure

   The .cif file can be converted by small programs to generate the fort.34 file, e.g. via:

   cif2crystal urea_113
   

   Note that the gas phase geometry will be automatically generated within a large unit cell (minimum distance to image of 12 Å) and the crystal geometry includes all symmetry transformations of the required space group.

   Hint

   Remember to activate the py27 environment to make cif2crystal work correctly.

2. Calculate the electronic energies at the TPSS-D3^ATM(BJ)/600 eV level. How large is the London dispersion contribution to the lattice energy?

   Technical procedure

   Fully optimized (TPSS-D3^ATM(BJ)) geometries for the urea molecule and the urea solid are available. Adjust the Brillouin sampling to a \(k\) grid of approximately 1/35 Bohr^-1. Use VASP 5.4 to calculate a TPSS-D3^ATM(BJ)/600 eV single-point energy for the gas and solid phase. The calculations can be done via the CRYSTAL 14 interface. Note that due to some internal restructuring for the VASP calculation, the interface broke. You can nevertheless obtain the correct energy using the command below to generate the vasp input file, do the VASP calculation by simply typing in vasp and then collect the results by again typing the command below.

   crystal_interface < INPUT > crystal.out
   

3. Calculate the free energy corrections at the semiempirical DFTB3-D3 level. Compare the electronic energy of the tight-binding model to the DFT one.

   Technical procedure

   1. Re-optimize the TPSS-D3^ATM(BJ) geometries of both phases with the CRYSTAL 14 program and the DFTB3-D3 Hamiltonian. Modify the INPUT file for the solid to sample the \(k\)-space identical to 1. Start the program by:

      crystal_interface < INPUT > crystal.out
      

   2. Calculate the phonon spectrum at the \(\Gamma\)-point with identical setup at the optimized geometry. In order to include low-lying (acoustic) modes of the solid, a supercell has to be generated (you find proper placeholders in the INPUT file). Increase the supercell until the energy contributions are converged within 0.1 kcal⋅mol^-1. The contributions can be taken from the output blocks in the calculation output entitled HARMONIC VIBRATIONAL CONTRIBUTIONS and THERMODYNAMIC FUNCTIONS.

4. Compare the calculated sublimation enthalpy with the experimental one. Explain why no thermal effects are included in the DFT energy and how this energy could be improved. Which significant approximations are applied in the overall model for the enthalpy correction?