In [ ]:
%pylab inline
import sys
from scipy.linalg import eigh
import moldesign as mdt
from moldesign import units as u
from moldesign.interfaces import nbo_interface
from itertools import product

In [ ]:
MOLNAME = 'conjugated_5'
QMMODEL = mdt.models.RHF(basis='6-31g')


First, we build a symmetric molecule.

In [ ]:
mol = mdt.from_smiles('C=CC#CC#CC#CC#CC=C')

In [ ]:

In [ ]:
mol.write(filename='dexterMols/%s.sdf' % MOLNAME) # TODO: shouldn't need keyword name
print 'wrote to ','dexterMols/%s.sdf'%MOLNAME


Next, we do a small - not necessarily converged - minimization to remove any energetic impossibilities.

TODO: use steepest descent for forces > 1 eV/Å; bfgs otherwise TODO: output tables not text TODO: visualize while running

In [ ]:
mol ='dexterMols/%s.sdf' % MOLNAME, format='sdf')

In [ ]:

In [ ]:
traj = mol.minimize(method='gradient descent', nsteps=40, frame_interval=1)
t2 = mol.minimize(method='bfgs', nsteps=60, frame_interval=2)  # TODO: monospace log fonts

In [ ]:

In [ ]:

In [ ]:
print 'Wrote dexterMols/%s_minimized.sdf' % MOLNAME
mol.write(filename='dexterMols/%s_minimized.sdf' % MOLNAME) # TODO: shouldn't need keyword "filename"

Setting up the data analysis

Next, we set up some involved analysis functions for the orbitals. This requires diving deeply into the results of the electronic structure.

In [ ]:
def get_valence_orbs(nboset):
    """Gives a truncated set of NBO orbitals with only bonding
    and anti-bonding orbitals on heavy atoms
    val_orbs = []
    for orb in nboset:
        if orb.nbotype in ('core','rydberg'): continue            
        isH = False
        for atom in orb.atoms:
            if atom.elem == 'H':
                isH = True
        if isH: continue
    valence = bb.orbitals.MolecularOrbitals(val_orbs,wfn=nboset.wfn,
    return valence

Triplet CIS (Configuration Interaction with Single excitations)

Now, we'll set up a hamiltonian in a basis of triplet-spin, singly-excited excitons. Each exciton is described in terms of the NLMO basis. Specifically, the basis function $|i,a>$ describes an electronic configuration in which spin-up electron was removed from excited NLMO orbital $i$, and a spin-down electron was created in NLMO anti-bonding orbital $a$.

The Hamiltonian in this basis is given by $$<i,a|H|j,b> = \delta_{ij} F_{ba} - \delta_{ab} F_{ij} - (ij|ab),$$ where $F_{xy}$ is the Fock matrix element between orbitals $x$ and $y$ and $(ij|ab)$ is the electron repulsion integral (ERI) between the four NLMO orbitals. We'll assemble an object to calculate these elements on the fly (since the matrix is really big and sparse, even for our limited set of orbitals).

In [ ]:
class TripletHamiltonian(object):
    def __init__(self, orbitals, eris=None, fock=None):
        self.orbs = orbitals
        if eris is None:
            print 'Computing 2e integrals ...',
            self.eri = orbitals.wfn.aobasis.get_eris_in_basis(self.orbs)
            print 'done.'
            self.eri = eris
        if fock is None:
            self.fock = orbitals.fock.defunits()
            self.fock = fock.defunits()
        = np.identity( len(orbitals) )
    def __getitem__(self,item):
        i,a,j,b = item
        return ([i,j]*self.fock[b,a] -[b,a]*self.fock[i,j]
                - self.eri[i,j,a,b])
    def fock_elem(self,i,a,j,b):
        return[i,j]*self.fock[b,a] -[b,a]*self.fock[i,j]
    def fock_exch_mod(self,i,a,j,b):
        # exch_delta is 1 if i==j AND/OR a==b
        exch_delta =[i,j] +[a,b] -[i,j]*[a,b]
        return ([i,j]*self.fock[b,a] 
                - exch_delta *self.eri[i,j,a,b])
class TripletCISWfn(object):
    def __init__(self, orbitals, nstates=10, _build=True, donor=None, acceptor=None):
        self.orbs = orbitals
        self.nstates = nstates

        self.occ = []
        self.virt = []

        for iorb,orb in enumerate(self.orbs):
            if orb.occupation == 2.0:
            elif orb.occupation == 0.0:
            else: raise ValueError('Orbital must have exactly 2 or 0 occupation')
        self.num_occ = len(self.occ)
        self.configs = [(i,j) for i,j in product(self.occ,self.virt)]
        self.config_index = {cfg:i for i,cfg in enumerate(self.configs)}
        self.nconfig = len(self.configs)
        self.config_orbs = [(self.orbs[i], self.orbs[j]) for i,j in self.configs]
        self.hamil = None
        if _build:
        if donor is not None and acceptor is not None:
            self.set_donor_acceptor(donor, acceptor)
        self.bonds, self.bond_dists = self.calc_bond_distances()
        self.bond_idxs = {b:i for i,b in enumerate(self.bonds)}
        self.bond_orbs = {}
        for iorb, orb in enumerate(self.orbs):
            if is not None:
    def build(self):
        if self.hamil is None: self.hamil = TripletHamiltonian(self.orbs)
        self.H, self.exch = self.buildh()
        self.evals = self.evecs = None
        self.state_coeffs, self.state_energies = self.diagonalize(self.nstates)
        self.dexter_coupling = (self.state_energies[1]-self.state_energies[0])/2.0
    def set_donor_acceptor(self, donor_tuple, acceptor_tuple):
        self.donor = tuple(donor_tuple)
        self.acceptor = tuple(acceptor_tuple)
        self.dstate = ciswfn.config_index[self.donor]
        self.astate = ciswfn.config_index[self.acceptor]
        self.etstate = ciswfn.config_index[self.donor[0],self.acceptor[1]]
        self.htstate = ciswfn.config_index[self.acceptor[0],self.donor[1]]
    def __getitem__(self, item):
        if len(item) == 2:
            return self.H[item]
            assert len(item) == 4
            return self.hamil[item]
    def buildh(self):
        """Use numpy indexing magic to make this super-fast"""
        i,a = zip(*self.configs)
        j,b = zip(*self.configs)
        ia = np.ix_(i,a)
        jb = np.ix_(j,b)
        H = self.hamil[ia[0],ia[1],jb[0].T,jb[1].T].defunits()
        exch = self.hamil.eri[ia[0],ia[1],jb[0].T,jb[1].T].defunits()
        return H, exch
    def shift_da_energies(self,shift):
        A positive shift will move the donor and acceptor orbitals away from the fermi level
        (i.e. moves virtual orbitals up by +shift, moves occupied orbitals down by -shift)
        newwfn = TripletCISWfn(self.orbs, nstates=self.nstates, _build=False,
                               donor=self.donor, acceptor=self.acceptor)
        newwfn.hamil = TripletHamiltonian(self.orbs,
                                          eris=self.hamil.eri, fock=self.hamil.fock.copy())
        for o in (self.donor, self.acceptor):
            newwfn.hamil.fock[o[0],o[0]] -= shift
            newwfn.hamil.fock[o[1],o[1]] += shift
        return newwfn
    def H_no_exchange(self):
        i,a = zip(*self.configs)
        j,b = zip(*self.configs)
        ia = np.ix_(i,a)
        jb = np.ix_(j,b)
        H = self.hamil.fock_elem(ia[0],ia[1],jb[0].T,jb[1].T).defunits()
        for idx in xrange(len(self.configs)): H[idx,idx] = self.H[idx,idx]
        return H
    def diagonalize(self, nstates):
        Solve for the lowest nstates states
        #print 'diagonalizing ...',;sys.stdout.flush()
        evals, evecs = eigh(self.H, eigvals=(0,nstates-1))
        self.evals = evals * self.H.get_units()
        self.evecs = evecs.T
        #print 'done.'
        return self.evecs, self.evals

    def kill_bonds(self, bondlist):
        killbonds = set(bondlist)
        mediating_configs = []
        for cfg,(orbh,orbe) in zip(self.configs,self.config_orbs):
            if cfg in (self.donor, self.acceptor): continue
            elif in killbonds: continue
            elif in killbonds: continue
            else: mediating_configs.append(cfg)
        bridge_ids = [self.config_index[x] for x in mediating_configs]
        return self.analyze(bridge=bridge_ids)
    def analyze(self, H=None, bridge=None):
        if H is None:
            H = self.H
        if bridge is None: 
            bridge = [i for i in xrange(self.nconfig) if i not in (self.dstate, self.astate)]
        newbridge = [b for b in bridge if b not in (self.etstate,self.htstate)]
        result = bb.utils.DotDict()           
        result.TET_G0 = self.lowdin_mediated_coupling(self.dstate, 
                                                      self.astate, bridge=bridge, H=H,
        result.HT_G0 = self.lowdin_mediated_coupling(self.dstate, 
                                                     bridge=newbridge, H=H,
                                                     direct=True, resonant=False)
        result.ET_G0 = self.lowdin_mediated_coupling(self.dstate, 
                                                     self.etstate, bridge=newbridge, H=H,
                                                     direct=True, resonant=False) = H[self.dstate,self.astate].defunits()
        result.CT = 2.0*(result.HT_G0 * result.ET_G0)/(
        result.scholes = result.CT +
        result.BE = self.lowdin_mediated_coupling(self.dstate, 
                                                  self.astate, bridge=newbridge, H=H,
                                                     direct=False, resonant=False)
        result.BE_subt = result.TET_G0 - result.scholes
        return result
    def lowdin_mediated_coupling(self, istate,fstate,H=None,
                                 bridge=None,E=None,direct=False, resonant=True):
        Calculates first-order superexchange coupling using first order perturbation theory
        :param istate: index of initial state
        :param fstate: index of the final state
        :param E: tunneling energy (default: H[iState,iState])
        :param direct: include the 0th-order contribution H[iState,fState]
        #Make sure calculation is reasonable
        if H is None: H = self.H
        if E is None: E=H[istate,istate]
        if resonant: assert np.allclose(H[istate,istate],H[fstate,fstate],rtol=1e-3)
        dimension = H.shape[0]
        assert H.shape[0] == H.shape[1]
        if bridge is None:
            bridge = [i for i in xrange(H.shape[0]) if i not in (istate,fstate)]
        assert istate not in bridge
        assert fstate not in bridge

        #Diagonalize subspace
        bridge_indices = np.array(bridge)
        bridge_selector = np.ix_(bridge_indices,bridge_indices)
        bridge_H = H[bridge_selector]
        ssevals,ssevecs = eigh(bridge_H)
        ssevals = ssevals * bridge_H.get_units()

        #Add it up
        lhs = H[istate,bridge_indices].dot(ssevecs)
        rhs = H[fstate,bridge_indices].dot(ssevecs)
        coup = (lhs/(E-ssevals)).dot(rhs)

        if direct: coup += H[istate,fstate]
        return coup
    def calc_bond_distances(self):
        bonds = set()
        for orb in self.orbs: bonds.add(
        bonds = list(bonds)
        nbonds = len(bonds)
        dists = np.inf*np.ones((nbonds,nbonds), dtype='int')
        for ibond, b in enumerate(bonds):
            for jbond in xrange(ibond+1,nbonds):
                nbr = bonds[jbond]
                if set([b.a1,b.a2]).intersection([nbr.a1,nbr.a2]):
                    dists[ibond, jbond] = dists[jbond, ibond] = 1
        for kbond in xrange(nbonds):
            for ibond in xrange(nbonds):
                for jbond in xrange(nbonds):
                    if dists[ibond,jbond] > dists[ibond, kbond] + dists[kbond, jbond]:
                        dists[ibond, jbond] = dists[jbond, ibond] = \
                        dists[ibond, kbond] + dists[kbond, jbond]
        return bonds, dists
    def only_nn_exchange(self, maxdist=1):
        H = self.H_no_exchange()
        occpairs = []
        virtpairs = []
        for iorb,orb1 in enumerate(self.orbs):
            if is None: continue
            for jorb in xrange(iorb,len(self.orbs)):
                orb2 = self.orbs[jorb]
                if is None: continue
                if orb1.occupation != orb2.occupation: continue
                ib = self.bond_idxs[]
                jb = self.bond_idxs[]
                if self.bond_dists[ib,jb] > maxdist: continue
                if orb1.occupation == 2.0:
                    assert orb1.occupation == 0.0
        for i,j in occpairs:
            for a,b in virtpairs:
                if (i==j) and (a==b): continue  # don't modify diagonals
                exch = self.hamil.eri[i,j,a,b]
                icfg = self.config_index[i,a]
                jcfg = self.config_index[j,b]
                H[icfg,jcfg] -= exch
                H[jcfg,icfg] -= exch
                if (i==j) or (a==b): continue
                icfg = self.config_index[i,b]
                jcfg = self.config_index[j,a]
                H[icfg,jcfg] -= exch
                H[jcfg,icfg] -= exch
        return H

Running the analysis

NLMO Orbital rotation

We'll rotate the molecular orbitals into a more useful basis - Natural Localized Molecular Orbitals, which reproduce Lewis-structure electron bonding.

In [ ]:
#!ls dexterMols/

In [ ]:

In [ ]:
mol ='dexterMols/%s_minimized.sdf' % MOLNAME)

In [ ]:

In [ ]:
nlmo_v = get_valence_orbs(mol.electronic_state.orbitals.NLMO)
mol.electronic_state.orbitals.nlmo_valence = nlmo_v

In [ ]:

Several sets of molecular orbitals should be available in the dropdown menu. We're interested specifically in the NLMOs, and even more specifically the truncated set that includes only heavy atom bonding/anti-bonding orbitals. It's important to verify that the localized donor and acceptor orbitals (the π and π* orbitals on the ethylenic moieties) are exactly degenerate with one another.

Angle scan

In [ ]:
rot_atoms = [mol.atoms[i] for i in (0,1,2,3)]
scan = bb.trajectory.Trajectory(mol)
mol.constrain_dihedral(*[mol.atoms[i] for i in (7,6,5,4)])

In [ ]:
wfns = []
mins = []
for deg in np.arange(-180,1.0,30)*u.degrees:
    print 'Dihedral:', deg
    bb.set_dihedral(*rot_atoms, theta=deg)
    mol.constraints[-1].value = deg
    scan.new_frame(annotation='dihedral:%s' % deg)
    nlmo_v = get_valence_orbs(mol.electronic_state.orbitals.NLMO)

In [ ]:

In [ ]:

Calcuating the overall coupling

Now, we build the CIS valence Hamiltonian

In [ ]:

In [ ]:

In [ ]:
ciswfn = TripletCISWfn(nlmo_v)
print ciswfn.dexter_coupling

Coupling pathway analysis

Select the donor and acceptor bonds

In [ ]:
da_selector = bb.ui.SelectionBuilder(mol)

In [ ]:
dbond, abond = da_selector.selected_bonds
d_idx = [None,None]
a_idx = [None, None]
for i,orb in enumerate(nlmo_v):
    if orb.nbotype != 'pi': continue
    if == dbond: arr = d_idx
    elif == abond: arr = a_idx
    else: continue # this isn't the pi-bond we're looking for
    if orb.occupation==2.0: arr[0] = i
    else: arr[1] = i
ciswfn.set_donor_acceptor(d_idx, a_idx)
tunneling_E = (ciswfn.state_energies[1]+ciswfn.state_energies[0])/2.0
TET_GE = ciswfn.lowdin_mediated_coupling(ciswfn.dstate, ciswfn.astate, direct=True, E=tunneling_E).defunits()
print 'Actual coupling', ciswfn.dexter_coupling
print 'Perturbative with full tunneling energy',TET_GE

In [ ]:

In [ ]:
class Slicer(list):
    def __getattr__(self, attr):
        return u.to_units_array([getattr(x,attr) for x in self])

Shift the D/A energies

These plots show how the coupling pathways are affected by the energetics of the donor and bridge. We apply an energy shift to the donor/acceptor orbitals - a positive shift moves them

In [ ]:
energies = np.arange(3,9,0.75) * u.eV
plot(ciswfn.hamil.fock.diagonal()[:ciswfn.num_occ], marker='_', ls='none',
     markersize=20.0, markeredgewidth=3.0, color='black')
plot(ciswfn.hamil.fock.diagonal()[ciswfn.num_occ:], marker='_', 
     ls='none', markersize=20.0, markeredgewidth=3.0, color='blue')

de = list()
results = Slicer()
no_exchange = Slicer()
for ee in energies:
    r = ciswfn.shift_da_energies(-ee)
    plot(r.hamil.fock.diagonal()[:ciswfn.num_occ], marker='_', ls='none',
     markersize=20.0, markeredgewidth=1.5, color='black')
    plot(r.hamil.fock.diagonal()[ciswfn.num_occ:], marker='_', 
     ls='none', markersize=20.0, markeredgewidth=1.5, color='blue')

grid(); xlabel('orbital index ->'); ylabel('orbital energy $F_{ii}$')
title('Shifted orbital energies')
#plot(energies, results.TET_G0, label='total (w/exchange)')
#plot(energies, no_exchange.TET_G0, label='total (no exchange)')
#grid(); xlabel('shift / eV'); ylabel('Coupling / eV')
#legend(); title(MOLNAME)
#plot(energies, results.scholes, label='scholes (w/exchange)', color='darkorange',lw=2)
#plot(energies, no_exchange.scholes, label='scholes (no exchange)', color='red')
#plot(energies, results.BE, label='BE (w/exchange)', color='darkblue', lw=2)
#plot(energies, no_exchange.BE, label='BE (no exchange)', color='blue')
#grid(); xlabel('shift / eV'); ylabel('Coupling / eV')
#legend(); title(MOLNAME+' without cross-exciton exchange')

plot(energies, results.TET_G0, label='total', lw=3, color='#0600ff')
plot(energies, results.BE, label='BE',lw=2, color='black')
plot(energies, results.scholes, label='scholes',lw=2, color='#00ff0b')
plot(energies,, label='$(da|d^*a^*)$', color='#ff0005', ls='--',lw=2)

legend(); title(MOLNAME)
grid(); xlabel('shift / eV'); ylabel('Coupling / eV')

In [ ]:
results = Slicer()
energies = np.arange(-2.5,5.0,1.0) * u.eV

nn_results = [Slicer() for i in xrange(4)]
for ee in energies:
    r = ciswfn.shift_da_energies(-ee)
    for i in xrange(4):
        H_nn = r.only_nn_exchange(maxdist=i)

In [ ]:
plot(energies, nn_results[0].TET_G0, label='$T_{ET}$ (no exch)', color=cm.inferno.colors[15], lw=2)
plot(energies, nn_results[1].TET_G0, label='$T_{ET}$ (NN exch)', color=cm.inferno.colors[1*256/4], lw=2)
plot(energies, nn_results[2].TET_G0, label='$T_{ET}$ (2nd NN exch)', color=cm.inferno.colors[2*256/4], lw=2)
plot(energies, nn_results[3].TET_G0, label='$T_{ET}$ (3rd NN exch)', color=cm.inferno.colors[3*256/4], lw=2)
plot(energies, results.TET_G0, label='$T_{ET}$ (full)', color=cm.inferno.colors[-15], lw=4)
grid(); xlabel('shift / eV'); ylabel('Coupling / eV')
legend(); title(MOLNAME+' dexter coupling')

plot(energies, nn_results[0].CT, label='CT (no exch)', color=cm.inferno.colors[15], lw=2)
plot(energies, nn_results[1].CT, label='CT (NN exch)', color=cm.inferno.colors[1*256/4], lw=2)
plot(energies, nn_results[2].CT, label='CT (2nd NN exch)', color=cm.inferno.colors[2*256/4], lw=2)
plot(energies, nn_results[3].CT, label='CT (3rd NN exch)', color=cm.inferno.colors[3*256/4], lw=2)
plot(energies, results.CT, label='CT (full)', color=cm.inferno.colors[-15], lw=4)
grid(); xlabel('shift / eV'); ylabel('Coupling / eV')
legend(); title(MOLNAME+' coupling through CT states')

plot(energies, nn_results[0].BE, label='BE (no exch)', color=cm.inferno.colors[15], lw=2)
plot(energies, nn_results[1].BE, label='BE (NN exch)', color=cm.inferno.colors[1*256/4], lw=2)
plot(energies, nn_results[2].BE, label='BE (2nd NN exch)', color=cm.inferno.colors[2*256/4], lw=2)
plot(energies, nn_results[3].BE, label='BE (3rd NN exch)', color=cm.inferno.colors[3*256/4], lw=2)
plot(energies, results.BE, label='BE (full)', color=cm.inferno.colors[-15], lw=4)
grid(); xlabel('shift / eV'); ylabel('Coupling / eV')
legend(); title(MOLNAME+' Bridge exciton coupling')

Select the bridge bonds to turn off

In [ ]:
kill_selector = bb.ui.SelectionBuilder(mol)

In [ ]:

In [ ]:
# Find bridge orbital indices

In [ ]: