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')
In [ ]:
mol = mdt.from_smiles('C=CC#CC#CC#CC#CC=C')
mdt.ui.GeometryBuilder(mol)
In [ ]:
bb.symmetry.Symmetrizer(mol)
In [ ]:
mol.write(filename='dexterMols/%s.sdf' % MOLNAME) # TODO: shouldn't need keyword name
print 'wrote to ','dexterMols/%s.sdf'%MOLNAME
In [ ]:
mol = bb.read('dexterMols/%s.sdf' % MOLNAME, format='sdf')
In [ ]:
mol.set_energy_model(QMMODEL)
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 [ ]:
t2.draw_orbitals()
In [ ]:
bb.symmetry.Symmetrizer(mol)
In [ ]:
print 'Wrote dexterMols/%s_minimized.sdf' % MOLNAME
mol.write(filename='dexterMols/%s_minimized.sdf' % MOLNAME) # TODO: shouldn't need keyword "filename"
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
break
if isH: continue
val_orbs.append(orb)
valence = bb.orbitals.MolecularOrbitals(val_orbs,wfn=nboset.wfn,
orbtype=nboset.orbtype+'_valence')
return valence
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 ...',
sys.stdout.flush()
self.eri = orbitals.wfn.aobasis.get_eris_in_basis(self.orbs)
print 'done.'
sys.stdout.flush()
else:
self.eri = eris
if fock is None:
self.fock = orbitals.fock.defunits()
else:
self.fock = fock.defunits()
self.delta = np.identity( len(orbitals) )
def __getitem__(self,item):
i,a,j,b = item
return (self.delta[i,j]*self.fock[b,a] - self.delta[b,a]*self.fock[i,j]
- self.eri[i,j,a,b])
def fock_elem(self,i,a,j,b):
return self.delta[i,j]*self.fock[b,a] - self.delta[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 = self.delta[i,j] + self.delta[a,b] - self.delta[i,j]*self.delta[a,b]
return (self.delta[i,j]*self.fock[b,a]
- self.delta[b,a]*self.fock[i,j]
- 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:
self.occ.append(iorb)
elif orb.occupation == 0.0:
self.virt.append(iorb)
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: self.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 orb.bond is not None:
self.bond_orbs.setdefault(orb.bond,[]).append(iorb)
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]
else:
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
newwfn.build()
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 orbh.bond in killbonds: continue
elif orbe.bond 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,
direct=True)
result.HT_G0 = self.lowdin_mediated_coupling(self.dstate,
self.htstate,
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)
result.direct = H[self.dstate,self.astate].defunits()
result.CT = 2.0*(result.HT_G0 * result.ET_G0)/(
H[self.dstate,self.dstate]-H[self.etstate,self.etstate])
result.scholes = result.CT + result.direct
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(orb.bond)
bonds = list(bonds)
nbonds = len(bonds)
dists = np.inf*np.ones((nbonds,nbonds), dtype='int')
np.fill_diagonal(dists,0)
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 orb1.bond is None: continue
for jorb in xrange(iorb,len(self.orbs)):
orb2 = self.orbs[jorb]
if orb2.bond is None: continue
if orb1.occupation != orb2.occupation: continue
ib = self.bond_idxs[orb1.bond]
jb = self.bond_idxs[orb2.bond]
if self.bond_dists[ib,jb] > maxdist: continue
if orb1.occupation == 2.0:
occpairs.append((iorb,jorb))
else:
assert orb1.occupation == 0.0
virtpairs.append((iorb,jorb))
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
In [ ]:
#!ls dexterMols/
In [ ]:
#MOLNAME =
In [ ]:
mol = bb.read('dexterMols/%s_minimized.sdf' % MOLNAME)
mol.draw()
In [ ]:
mol.set_energy_model(QMMODEL)
mol.calculate()
mol.electronic_state.run_nbo()
In [ ]:
nlmo_v = get_valence_orbs(mol.electronic_state.orbitals.NLMO)
mol.electronic_state.orbitals.nlmo_valence = nlmo_v
In [ ]:
mol.draw_orbitals()
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.
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)])
mol.constrain_dihedral(*rot_atoms)
In [ ]:
mol.set_energy_model(bb.models.PySCFPotential(theory='rhf',basis='sto-3g'))
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
mins.append(mol.minimize(nsteps=30,frame_interval=1))
mol.electronic_state.run_nbo()
scan.new_frame(annotation='dihedral:%s' % deg)
nlmo_v = get_valence_orbs(mol.electronic_state.orbitals.NLMO)
wfns.append(TripletCISWfn(nlmo_v))
In [ ]:
In [ ]:
mins[0].draw3d()
In [ ]:
bb.ui.GeometryBuilder(mol)
In [ ]:
bb.geometry.dihedral_gradient(*rot_atoms)
In [ ]:
ciswfn = TripletCISWfn(nlmo_v)
print ciswfn.dexter_coupling
In [ ]:
da_selector = bb.ui.SelectionBuilder(mol)
da_selector
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 orb.bond == dbond: arr = d_idx
elif orb.bond == 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 [ ]:
ciswfn.analyze()
In [ ]:
class Slicer(list):
def __getattr__(self, attr):
return u.to_units_array([getattr(x,attr) for x in self])
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')
de.append(r.H[ciswfn.dstate,ciswfn.dstate])
no_exchange.append(r.analyze(r.H_no_exchange()))
results.append(r.analyze())
grid(); xlabel('orbital index ->'); ylabel('orbital energy $F_{ii}$')
xlim(-0.5,ciswfn.num_occ+1)
title('Shifted orbital energies')
figure()
#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)
#figure()
#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')
figure()
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, results.direct, 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)
nn_results[i].append(r.analyze(H_nn))
results.append(r.analyze())
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')
figure()
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')
figure()
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')
In [ ]:
kill_selector = bb.ui.SelectionBuilder(mol)
kill_selector
In [ ]:
ciswfn.analyze()
In [ ]:
# Find bridge orbital indices
ciswfn.kill_bonds(kill_selector.selected_bonds)
In [ ]: