charge transfer
Contributed by Ramon Crehuet for pDynamo version 1.7.2.
The code on this page is available here.
Introduction
The Marcus theory allows the evaluation of the rate of charge transfer (CT) between two molecules or two fragments of the same molecule, under certain circumstances (the so-called temperature activated, non-adiabatic limit). This theory and its subsequent generalizations have found several applications, both in biological and in materials science contexts. One important ingredient of the CT rate expression is the so-called "transfer integral" ( tij ) between orbital i on one molecule and orbital j on the other molecule. Typically, these would be the two HOMOs (for transfer of positively charged "holes") or the two LUMOs (for transfer of negatively charge electrons). The "reorganization energy" and the "site energy difference" between the two molecules are two other important ingredients, which however will not concern us here.
The aim of this section is to show how to compute the CT transfer integrals with pDynamo, using a pair of ethylene molecules as a simple example. Among many possible strategies, I have chosen the dimer projection (DIPRO) method described by Baumeier, Kirkpatrick and Andrienko (PCCP 2010, 12, 11103), which is approximate but has the advantage of being quite simple. A very similar method had also been described earlier, by Bredas and coworkers (JACS 2006, 128, 9882). It can be applied to both semiempirical and ab initio Hamiltonians (with some restrictions on the latter, due to limitations of the DFT code within pDynamo 1.7.2 which should be remedied soon). Time and space limitations do not allow a discussion of this approach. Those interested in the details should consult the original BKA paper. Here I use the "no counterpoise" version of the DIPRO method, whereby the input to the calculation is represented by the molecular orbitals of the two neutral, isolated molecules and of their supramolecular complex.
The code examples in this section also demonstrate how to extract and manipulate some pDynamo arrays, which in this case contain the AO overlap matrix and the molecular orbital coefficients (two matrices or 2D arrays) and the orbitals energies (a vector or 1D array) of a QC system. Note that the AO overalp matrix is unity in the MNDO approximation, whereas it is a general symmetric matrix in the ab initio case.
Extracting orbitals and orbital energies
The code below should be almost self-explanatory:
from pCore import Clone
from pBabel import XYZFile_ToSystem
from pMolecule import ElectronicState, QCModelMNDO, QCModelDFT
from pMoleculeScripts import ConjugateGradientMinimize_SystemGeometry
doOptimize = True
qcmodel = [ QCModelMNDO( "am1", isSpinRestricted = True, keepOrbitalData = True ) , \
QCModelDFT( functional = "blyp", orbitalBasis = "631gs", isSpinRestricted = True, keepOrbitalData = True ) ]
for qc in qcmodel:
mol1 = XYZFile_ToSystem( "Ethylene1.xyz" )
mol1.DefineQCModel( qc )
mol1.Summary ( "Summary for molecule 1" )
energy1 = mol1.Energy ( doGradients = True )
if doOptimize:
xyz1 = Clone( mol1.coordinates3 )
ConjugateGradientMinimize_SystemGeometry( mol1, logFrequency=10, maximumIterations=200, rmsGradientTolerance=1.0 )
mol1.coordinates3.Superimpose( xyz1 )
rms1 = xyz1.RMSDeviation( mol1.coordinates3 )
energy1 = mol1.Energy()
(eps1, homo1, lumo1) = mol1.energyModel.qcModel.OrbitalEnergies( mol1.configuration )
eps1.Print( title = "Orbital Energies" )
mos1 = mol1.energyModel.qcModel.Orbitals( mol1.configuration )
mos1.Print( title = "Orbitals" )
Note: using pDynamo-1.7.2, this program runs smoothly for the semiempirical part, but it fails in the DFT part. The reason is that extracting DFT orbitals had not yet been implemented, yet. I have been able to run the code thanks to a patch, which should become available to everyone with the next release of pDynamo.
Charge transfer calculation: driver program
from pCore import Clone, YAMLPickle, YAMLUnpickle, Real2DArray, Real1DArray, SymmetricMatrix
from pBabel import XYZFile_ToSystem
from pMolecule import QCModelMNDO, QCModelDFT
from pMoleculeScripts import MergeByAtom, ConjugateGradientMinimize_SystemGeometry
from TM import TransferMatrix
rhf = QCModelMNDO ( "rm1", isSpinRestricted = True, keepOrbitalData=True )
#dft = QCModelDFT ( functional = "lda", orbitalBasis = "sto3g", isSpinRestricted = True, keepOrbitalData=True )
# . Molecule 1
mol1 = XYZFile_ToSystem( "Ethylene1.xyz" )
mol1.DefineQCModel( rhf )
energy1 = mol1.Energy ( doGradients = True )
YAMLPickle( "Ethylene1.yaml", mol1 )
(eps1, homo1, lumo1) = mol1.energyModel.qcModel.OrbitalEnergies( mol1.configuration )
print "HOMO and LUMO indexes for 1:", homo1, lumo1
eps1.Print( title = "Orbital Energies for 1" )
mos1 = mol1.energyModel.qcModel.Orbitals( mol1.configuration )
mos1[ : , homo1 ].Print( title = "HOMO Orbital for 1" )
mos1[ : , lumo1 ].Print( title = "LUMO Orbital for 1" )
# . Molecule 2
mol2 = XYZFile_ToSystem( "Ethylene2.xyz" )
mol2.DefineQCModel( rhf )
energy2 = mol2.Energy ( doGradients = True )
YAMLPickle( "Ethylene2.yaml", mol2 )
(eps2, homo2, lumo2) = mol2.energyModel.qcModel.OrbitalEnergies( mol2.configuration )
print "HOMO and LUMO indexes for 2:", homo2, lumo2
eps2.Print( title = "Orbital Energies for 2" )
mos2 = mol2.energyModel.qcModel.Orbitals( mol2.configuration )
mos2[ : , homo2 ].Print( title = "HOMO Orbital for 2" )
mos2[ : , lumo2 ].Print( title = "LUMO Orbital for 2" )
# . Dimer (molecule 1 + molecule 2)
dimer = MergeByAtom( mol1, mol2 )
dimer.DefineQCModel( rhf )
energy12 = dimer.Energy ( doGradients = True )
YAMLPickle( "EthyleneDimer.yaml", dimer )
(eps12, homo12, lumo12) = dimer.energyModel.qcModel.OrbitalEnergies( dimer.configuration )
print "HOMO and LUMO indexes for dimer:", homo12, lumo12
eps12.Print( title = "Orbital Energies for dimer" )
mos12 = dimer.energyModel.qcModel.Orbitals( dimer.configuration )
interEnergy = energy12 - energy1 - energy2
print "Interaction Energy=", interEnergy
# . The calculation of CT Hamiltonian involves a symmetric (Lowdin) orthogonalization of the overlap matrix between
# . orbitals of the two fragments.
S = Clone ( dimer.configuration.qcState.overlap )
# . The real calculation is carried out by the TransferMatrix function
( Hh, Hl ) = TransferMatrix(mos1, homo1, lumo1, mos2, homo2, lumo2, mos12, eps12, S)
Hh.Print ( title = "HOMO-HOMO effective Hamiltonian for hole+ transport" )
Hl.Print ( title = "LUMO-LUMO effective Hamiltonian for elec- transport" )
Charge transfer calculation: the TransferMatrix function
from pCore import Real1DArray, Real2DArray, SymmetricMatrix
from pMolecule import UNITS_ENERGY_HARTREES_TO_ELECTRON_VOLTS
def TransferMatrix(mos1, homo1, lumo1, mos2, homo2, lumo2, mos12, eps12, S):
"""Returns the effective Hamiltonian matrices for hole (homo-homo) and electron (lumo-lumo) transfer
Energies are in eV.
"""
if (not mos1.isSquare or not mos2.isSquare or not mos12.isSquare):
return None
n1 = mos1.rows
n2 = mos2.rows
n12 = mos12.rows
if ( n1+n2 != n12 ):
return None
vec1 = Real2DArray(n12,n12)
vec2 = Real2DArray(n12,n12)
mos1.CopyTo( vec1[ :n1 , :n1 ] )
mos2.CopyTo( vec2[ n1: , n1: ] )
homo2 += n1
lumo2 += n1
vec1[ : , homo1 ].Print( title = "0-padded HOMO orbital for 1" )
vec1[ : , lumo1 ].Print( title = "0-padded LUMO orbital for 1" )
vec2[ : , homo2 ].Print( title = "0-padded HOMO orbital for 2" )
vec2[ : , lumo2 ].Print( title = "0-padded LUMO orbital for 2" )
temp = Real1DArray(n12)
gamma1h = Real1DArray(n12)
gamma1l = Real1DArray(n12)
gamma2h = Real1DArray(n12)
gamma2l = Real1DArray(n12)
# . Evaluate the HOMO1-HOMO2 and LUMO1-LUMO2 overlaps and the gamma vectors (for HOMO1, HOMO2, LUMO1, LUMO2)
if (S is None):
# . MNDO case: the AO overlap matrix S is the identity
print "S is None (MNDO case)"
Shh = vec1[ : , homo1 ].Dot( vec2[ : , homo2 ] )
Sll = vec1[ : , lumo1 ].Dot( vec2[ : , lumo2 ] )
print "HOMO1-HOMO2 Overlap = ", Shh
print "LUMO2-LUMO2 Overlap = ", Sll
mos12.VectorMultiply( vec1[ : , homo1 ], gamma1h, transpose = True )
mos12.VectorMultiply( vec1[ : , lumo1 ], gamma1l, transpose = True )
mos12.VectorMultiply( vec2[ : , homo2 ], gamma2h, transpose = True )
mos12.VectorMultiply( vec2[ : , lumo2 ], gamma2l, transpose = True )
else:
# . Ab initio DFT case: the AO overlap matrix S is not the identity
print "S isn't None (DFT case)"
S.VectorMultiply( vec2[ : , homo2 ] , temp )
Shh = vec1[ : , homo1 ].Dot( temp )
S.VectorMultiply( vec2[ : , lumo2 ] , temp )
Sll = vec1[ : , lumo1 ].Dot( temp )
print "HOMO1-HOMO2 Overlap = ", Shh
print "LUMO2-LUMO2 Overlap = ", Sll
S.VectorMultiply( vec1[ : , homo1 ] , temp )
mos12.VectorMultiply( temp, gamma1h, transpose = True )
S.VectorMultiply( vec1[ : , lumo1 ] , temp )
mos12.VectorMultiply( temp, gamma1l, transpose = True )
S.VectorMultiply( vec2[ : , homo2 ] , temp )
mos12.VectorMultiply( temp, gamma2h, transpose = True )
S.VectorMultiply( vec2[ : , lumo2 ] , temp )
mos12.VectorMultiply( temp, gamma2l, transpose = True )
# . HOMO1-HOMO2 matrix elements (before orthogonalization)
gamma1h.CopyTo( temp )
temp.Multiply( eps12 )
Eh11 = gamma1h.Dot( temp )
Eh12 = gamma2h.Dot( temp )
gamma2h.CopyTo( temp )
temp.Multiply( eps12 )
Eh22 = gamma2h.Dot( temp )
# . HOMO1-HOMO2 effective Hamiltonian matrix (after orthogonalization)
Hh = SymmetricMatrix(2)
norm = 1.0 - Shh**2
Hh[0,0] = 0.5 * ( Eh11+Eh22 - 2*Eh12*Shh + (Eh11-Eh22)*norm**0.5 ) / norm
Hh[1,1] = 0.5 * ( Eh11+Eh22 - 2*Eh12*Shh - (Eh11-Eh22)*norm**0.5 ) / norm
Hh[0,1] = ( Eh12 - 0.5*(Eh11+Eh22)*Shh ) / norm
Hh.Scale( UNITS_ENERGY_HARTREES_TO_ELECTRON_VOLTS )
# . LUMO1-LUMO2 matrix elements (before orthogonalization)
gamma1l.CopyTo( temp )
temp.Multiply( eps12 )
El11 = gamma1l.Dot( temp )
El12 = gamma2l.Dot( temp )
gamma2l.CopyTo( temp )
temp.Multiply( eps12 )
El22 = gamma2l.Dot( temp )
# . LUMO1-LUMO1 effective Hamiltonian matrix (after orthogonalization)
Hl = SymmetricMatrix(2)
norm = 1.0 - Sll**2
Hl[0,0] = 0.5 * ( El11+El22 - 2*El12*Sll + (El11-El22)*norm**0.5 ) / norm
Hl[1,1] = 0.5 * ( El11+El22 - 2*El12*Sll - (El11-El22)*norm**0.5 ) / norm
Hl[0,1] = ( El12 - 0.5*(El11+El22)*Sll ) / norm
Hl.Scale( UNITS_ENERGY_HARTREES_TO_ELECTRON_VOLTS )
# . Finished
return ( Hh, Hl )
This is the part of the code which really implements the DIPRO method by BKA. As can be seen, it can be cast in terms of a few matrix-matrix and matrix-vector multiplications.