About Issues Tutorials Documentation </span>

Example 2: Using MD sampling to calculate UV-Vis spectra

This notebook uses basic quantum chemical calculations to calculate the absorption spectra of a small molecule.

Author: Aaron Virshup, Autodesk Research
Created on: September 23, 2016
Tags: excited states, CASSCF, absorption, sampling



In [ ]:

    
%matplotlib inline
import numpy as np
from matplotlib.pylab import *

try: import seaborn  #optional, makes plots look nicer
except ImportError: pass

import moldesign as mdt
from moldesign import units as u

Single point
Sampling
Post-processing
Create spectrum

Single point

Let's start with calculating the vertical excitation energy and oscillator strengths at the ground state minimum (aka Franck-Condon) geometry.

Note that the active space and number of included states here is system-specific.



In [ ]:

    
qmmol = mdt.from_name('benzene')
qmmol.set_energy_model(mdt.models.CASSCF, active_electrons=6,
                       active_orbitals=6, state_average=6, basis='sto-3g')



In [ ]:

    
properties = qmmol.calculate()

This cell print a summary of the possible transitions.

Note: you can convert excitation energies directly to nanometers using Pint by calling energy.to('nm', 'spectroscopy').



In [ ]:

    
for fstate in range(1, len(qmmol.properties.state_energies)):
    excitation_energy = properties.state_energies[fstate] - properties.state_energies[0]
    
    print('--- Transition from S0 to S%d ---' % fstate )   
    print('Excitation wavelength: %s' % excitation_energy.to('nm', 'spectroscopy'))
    print('Oscillator strength: %s' % qmmol.properties.oscillator_strengths[0,fstate])

Sampling

Of course, molecular spectra aren't just a set of discrete lines - they're broadened by several mechanisms. We'll treat vibrations here by sampling the molecule's motion on the ground state at 300 Kelvin.

To do this, we'll sample its geometries as it moves on the ground state by:

Create a copy of the molecule
Assign a forcefield (GAFF2/AM1-BCC)
Run dynamics for 5 ps, taking a snapshot every 250 fs, for a total of 20 separate geometries.



In [ ]:

    
mdmol = mdt.Molecule(qmmol)
mdmol.set_energy_model(mdt.models.GaffSmallMolecule)
mdmol.minimize()



In [ ]:

    
mdmol.set_integrator(mdt.integrators.OpenMMLangevin, frame_interval=250*u.fs,
                     timestep=0.5*u.fs, constrain_hbonds=False, remove_rotation=True,
                     remove_translation=True, constrain_water=False)
mdtraj = mdmol.run(5.0 * u.ps)

Post-processing

Next, we calculate the spectrum at each sampled geometry. Depending on your computer speed and if PySCF is installed locally, this may take up to several minutes to run.



In [ ]:

    
post_traj = mdt.Trajectory(qmmol)
for frame in mdtraj:
    qmmol.positions = frame.positions
    qmmol.calculate()
    post_traj.new_frame()

This cell plots the results - wavelength vs. oscillator strength at each geometry for each transition:



In [ ]:

    
wavelengths_to_state = []
oscillators_to_state = []

for i in range(1, len(qmmol.properties.state_energies)):
    wavelengths_to_state.append(
        (post_traj.state_energies[:,i] - post_traj.potential_energy).to('nm', 'spectroscopy'))
    oscillators_to_state.append([o[0,i] for o in post_traj.oscillator_strengths])
    
for istate, (w,o) in enumerate(zip(wavelengths_to_state, oscillators_to_state)):
    plot(w,o, label='S0 -> S%d'%(istate+1),
         marker='o', linestyle='none')
xlabel('wavelength / nm'); ylabel('oscillator strength'); legend()

Create spectrum

We're finally ready to calculate a spectrum - we'll create a histogram of all calculated transition wavelengths over all states, weighted by the oscillator strengths.



In [ ]:

    
from itertools import chain
all_wavelengths = u.array(list(chain(*wavelengths_to_state)))
all_oscs = u.array(list(chain(*oscillators_to_state)))
hist(all_wavelengths, weights=all_oscs, bins=50)
xlabel('wavelength / nm')



In [ ]:

Example 2: Using MD sampling to calculate UV-Vis spectra

Contents

Single point

Sampling

Post-processing

Create spectrum