718
J. Chem. Theory Comput. 2010, 6, 718–726
Blending Determinism with Evolutionary Computing: Applications to the Calculation of the Molecular Electronic Structure of Polythiophene Kanchan Sarkar, Rahul Sharma, and S. P. Bhattacharyya* Department of Physical Chemistry, Indian Association for the CultiVation of Science, JadaVpur, Kolkata-700032, India Received October 13, 2009
Abstract: A density matrix based soft-computing solution to the quantum mechanical problem of computing the molecular electronic structure of fairly long polythiophene (PT) chains is proposed. The soft-computing solution is based on a “random mutation hill climbing” scheme which is modified by blending it with a deterministic method based on a trial single-particle density matrix [P(0)(R)] for the guessed structural parameters (R), which is allowed to evolve under a unitary transformation generated by the Hamiltonian H(R). The Hamiltonian itself changes as the geometrical parameters (R) defining the polythiophene chain undergo mutation. The scale (λ) of the transformation is optimized by making the energy [E(λ)] stationary with respect to λ. The robustness and the performance levels of variants of the algorithm are analyzed and compared with those of other derivative free methods. The method is further tested successfully with optimization of the geometry of bipolaron-doped long PT chains.
1. Introduction Exact solutions of the Schrodinger equation for manyelectron systems are impossible to obtain analytically. The solutions are therefore obtained by various approximation methods, which are still rather complicated for applications to large systems. With rapid and spectacular advances in digital computing, computational methods of electronic structure calculations have attracted serious attention in recent years.1-3 These methods have become practically essential for understanding large systems at the microscopic level. Finding efficient algorithms for handling large quantum systems even at an approximate level has become a challenging issue. Traditional electronic structure algorithms calculate eigenstates associated with discrete energy levels, and this leads to a diagonalization problem of the Hamiltonian matrix. The computational effort in traditional diagonalization methods scales as N3, where N is the dimension of the basis space. Usually one is interested in finding the lowest energy structure so that one has to simultaneously search the potential energy surface E(R) for locating the * Corresponding author fax: (91)332473 2805; e-mail: pcspb@ iacs.res.in.
global minimum while diagonalizing the Hamiltonian matrix H(R) at different geometries (R). In the present paper we have proposed a novel nondeterministic algorithm for locating the minimum energy structures of neutral or doped PT chains of 100, 150, and 200 thiophene rings. To be specific, we have proposed a variant of the “directed random mutation hill climbing (DRMHC) method”,14,15 which works with a randomly generated string of all the geometrical parameters (nuclear position variables, R) required to compute the energy and therefore the fitness of a neutral or doped PT molecule within the framework of a modified SSH effective π-electron Hamiltonian model.12,13 The geometry string {Ri} is allowed to undergo directed random mutations. The string of mutated geometry variables {Ri′} is used to define the Hamiltonian H({Ri′}), which acts as the generator of a unitary transformation, U(λ,{Ri′}). U(λ,{R′}) transforms a trial one-electron density matrix (P(0)) i into a “mutated” one-electron density matrix (P(0)′), which is used to compute the “energy” of the mutated structure and hence the fitness of the structure coded by the mutated geometry string. The parameter λ fixes the scale of the transformation which is optimized at each random mutation hill climbing step by making the energy stationary with
10.1021/ct900540d 2010 American Chemical Society Published on Web 02/01/2010
Blending Determinism with Evolutionary Computing
Figure 1. Fragment of the thiophene chain: (a) C-C bridging bond, (b) C-S single bond, (c) C-C double bond, (d) C-C single bond. Table 1. Parameters in the SSH Hamiltonian (for Polythiophene) Used in These Calculations param
value
description
A B D 0 R(C-C) 0 R(C-S) 0 R(C-C)B
123.6 eV 0.3776 au 7.814 au-1 1.557 au 1.782 au 1.557 au
Hamiltonian parameter Hamiltonian parameter Hamiltonian parameter C-C single bond length C-S bond length C-C bridging bond length
respect to λ. We note here that there have been a number of earlier attempts to calculate the one- and two particle-density matrices using similar unitary transformations8-10 on a trial density. We have, to the best of our knowledge, tried for the first time to blend it with an otherwise stochastic approach for computing the one-electron density matrix. We have compared the λ-optimized DRMHC algorithm with several alternative soft-computing methods of solving the problem. The algorithm is elaborately tested on long polythiophene chains, undoped as well as doped, and the results are briefly analyzed. The outline of the paper is as follows. In section 2 we present the algorithm in detail, while the results of the applications are described in section 3. Section 4 presents the concluding remarks.
2. Method Polythiophene (PT) oligomers have attracted the widespread attention of experimentalists and theoreticians because these molecules and their derivatives are capable of displaying metallic conduction under electron or hole doping.4-7,11-13 Neutral PTs are chemically stable, can be synthesized easily, and can be doped with dopants such as ClO4 and AsF5. The stability remains intact even after doping. The UV photoelectron spectrum5 of neutral PT shows the existence of two π-bands in the system, and the Fermi level is ∼1.2 eV above the valence band (VB). There are two nondegenerate classical resonating forms of PT in which the 2pz orbitals of the carbon atoms and the 3pz orbitals of the sulfur atoms interact to form a π-band, half of which is occupied by the electrons. Had the two forms been isoenergic, a solitonic mode of conduction would have been viable as the solitons have a tendency to separate because the structure between the defects is isoenergic with the structure outside. However, in PT, the degeneracy is weakly lifted, and as a consequence the solitonic mode of conduction is ruled out. It turns out that polarons and bipolarons are the most important excitations and charge-storage configurations in doped PT.6 While PT is a semiconductor with a band gap of about 2 eV, hole doping reduces the gap, and at high doping levels
J. Chem. Theory Comput., Vol. 6, No. 3, 2010 719
PT begins to show metallic properties. It is important therefore to understand how the electronic structure of neutral PT evolves under doping. It is expedient in this context to have algorithms that can locate the global minimum energy structures on the potential energy surface of PT oligomers, undoped or doped. There are several deterministic methodologies, mostly based on a gradient search, which have been explored for locating the minimum energy structures of PTs.6,7 Rarely, however, nondeterministic techniques have been explored in this context.11 Our primary concern here is to introduce such a hybrid technique of geometry optimization and calculation of electronic structure and test its workability with undoped and doped PTs as examples. PT is treated as a conjugated polyene and described by a tight binding model of the π-electronic system that includes only the nearest neighbor hopping interactions. The lattice dynamics is treated classically. This leads to a modified SSH effective π-electronic Hamiltonian H where12,13 H ) Hπe + Hlσ
(1)
The π-electronic Hamiltonian (Hπe ) is expressed as Hπe )
∑ Rna†nan + 21 ∑ Vnma†nam + hc n
(2)
n,m
and the lattice Hamiltonian (the σ-electronic framework Hamiltonian) is defined as Hlσ
)
∑ n
pn2 1 + 2Mn 2
∑ f(Rnm)
(3)
n,m
Rn is the self-energy of the carbon 2pz electron (3pz electron of sulfur) at the nth site, and Vnm represents the hopping interaction between the nth and mth sites. Vnm is parametrized in the form7,11 Vnm ) -Ae-Rnm/B ) 0
(|n - m| ) 1) (otherwise)
(4)
with two parameters, A and B. Rnm is the distance separating the nth and mth adjacent sites. We use frozen core approximation so that the kinetic energy of the lattice is zero, i.e. Hlσ )
1 2
∑ f(Rnm)
(5)
n,m
where 0 f(Rnm) ) -ADnm(Rnm - Rnm + B)e-Rnm/B 0 A, B, and Dnm are parameters of the model and Rnm is the 7,11 standard length of the n-m bond. We start by guessing randomly the M number of bond lengths required to describe the Nr-ring chain, forming and diagonalizing the N-dimensional (N ) 5Nr) π-electron Hamiltonian Hπe in the N-dimensional basis of the carbon 2pz and sulfur 3pz atomic orbitals of all the carbon and sulfur atoms of the chain. The diagonalization generates a set of N π molecular orbitals {φi}, Noc of which are occupied by the
720
J. Chem. Theory Comput., Vol. 6, No. 3, 2010
Sarkar et al.
π-electrons in pairs (Noc ) 3Nr). In general, we can write the π MOs {φi} as linear combinations of the atomic basis orbitals (χp), where N
φi ) Hπe φi
)
∑ nkεπk ) 2 ∑ Rnq0n
)
k)1
+2
n)1
∑ cp χp
p)1 εiπφi,
Noc
EGπ
(6)
i
i ) 1, 2, ..., N
επi represent the binding energies of the π molecular electrons. The linear expansion coefficients {cpi} define the elements of the charge density bond order matrix P through the relation
0 ∑ ∑ VnmPnm
n
