Microsoft Word - European_Journal_of_Chemistry_1_2_2010_142_161_119 European Journal of Chemistry 1 (2) (2010) 142‐161    European Journal of Chemistry  ISSN 2153‐2249 (Print) / ISSN 2153‐2257 (Online)  2010 EURJCHEM  DOI:10.5155/eurjchem.1.2.142‐161.119        European Journal of Chemistry  Journal homepage: www.eurjchem.com          Theory and computations of two‐photon absorbing photochromic  chromophores  Artëm E. Masunova,b,c,d,* and Ivan A. Mikhailova  a NanoScience Technology Center, University of Central Florida, Orlando, FL 32826, USA  b Department of Chemistry, University of Central Florida, Orlando, FL 32816, USA  c Department of Physics, University of Central Florida, Orlando, FL 32816, USA  d Florida Solar Energy Center, University of Central Florida, Cocoa, FL 32922, USA  *Corresponding author at: NanoScience Technology Center, University of Central Florida, Orlando, FL 32826, USA. Tel.: +1.321.6622724; fax: +1.407.8822819.  E­mail address: amasunov@gmail.com (A.E. Masunov).                REVIEW INFORMATION    ABSTRACT Received: 11 March 2010  Received in revised form: 20 May 2010  Accepted: 23 May 2010  Online: 30 June 2010  KEYWORDS    Exponential  growth  in  information  technology  generates  ever  increasing  amounts  of  data,  making recording density of  the storage media crucially  important. Two‐photon absorption  was proposed as a basis for high‐density multi‐layer technology for optical memory and logic  devices.  This  technology  suggests  to  use  polymers,  doped  with  photochromic  compounds  that undergo a reversible photoinduced isomerization, or photoswitching. In this review we  consider  recent  theoretical  works  and  benchmarking  studies  of  the  DFT‐based  methods,  capable to predict two‐photon absorption (2PA) and photochemical activity, Next we review  the applications of these methods to design a prototype molecule that combines the photon‐ mode  recording  property  of  photochromic  compounds  with  large  2PA  cross‐section.  We  conclude that a posteriori Tamm‐Dancoff approximation to the second order CEO approach in  Density  Functional  Theory  is  the  powerful  tool  for  both  quantitative  predictions  and  qualitative understanding of the excited state processes in photophysics and photochemistry.  We  also  emphasize  general  principles  for  the  rational  design  of  a  two‐photon  operated  photoswitch.  TDDFT  Organic photochemistry  Photocyclization  Photoinduced reactions   Two‐photon absorption  State‐to‐state transition dipoles    1. Introduction    Exponential growth in the information technology industry  generates  ever  increasing  amounts  of  data.  This  made  the  recording  density  of  storage  media  crucially  important,  and  prompted development of optical solutions as an alternative to  the  traditional magnetic  storage  devices.  Optical  data  storage  refers to systems that use laser light for recording information  as well as for readout. Most optical memory systems available  on  the  market  currently  utilize  inorganic  materials  and  are  based  on  magneto‐optical  [1,2]  and  phase‐change  [3]  effects.  Both effects are utilized in the heat‐mode recording, where the  light energy is converted into heat. This induces a magnetic or  structural phase transition, and changes physical properties of  the medium.  Until recently, organic materials have not been considered  as  viable  candidates  because  of  insufficient  reliability.  However,  the situation has changed after CD‐R (compact disk‐ recordable)  technology  was  introduced.  This  technology  is  based  on  organic  dyes  as  the  memory  medium,  and  uses  polymers,  doped  with  photochromophores  (the  molecules,  which  undergo  a  reversible  photoinduced  isomerization  or  photochromism)  [4].  Photochromism  is  a  non‐destructive  process  involving  a  light  initiated  rearrangement  of  chemical  bonds accompanied by change in color and other properties. An  example  of  photochromophores  is  diarylethene  compounds,  shown  on  Scheme  1.  They  undergo  photochromic  processes,  involving the photoinduced conrotatory opening and closing of  the  central bond,  and  fulfill  requirements needed  for material  applications, namely,  thermal stability and resistance to  linear  optical photofatigue [5].        Scheme 1. Diarylethene photoswitching X = O or X=S. h stands for photon  absorption,  – thermal process.   Photon‐mode  data  recording  is  expected  to  win  the  competition  with  magnetic  and  heat‐mode  optical  recording.  Photochromic  materials  are  now  considered  to  be  promising  candidates  for erasable memory media of  the next generation.  Another  important  property  of  photochromic  materials  is  an  ultrafast  change  of  the  state  upon  irradiation.  Photochromic  compounds  with  ultrafast  switching  capabilities  can  be  also  useful  for  optical  communication,  signal  processing,  and  interconnection  as  components  of  various  photonic  devices,  such as optical switches, variable frequency filters, attenuators,  and phase shifters.  Besides reversibility and ultrafast switching, photochromic  optical  data  storage  devices  (as  opposed  to  heat‐mode  ones)  possess  another  important  advantage.  They offer  a  possibility  of data recording using individually addressable 3‐dimensional  arrays of bits. Three‐dimensional design opens an opportunity  of  using  hundreds  of  layers  beneath  the  optical  disk  surface,  thus  increasing  storage capacity  two orders of magnitude and  more.  In  order  to  increase  a  storage  density  one  needs  to  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  143  reduce  the parasite  cross‐talk between  these  recorded  layers;  to  this  end  Rentzepis  et  al.  suggested  using  two‐photon  absorption of light [6].                Figure 1.  One‐photon  (blue  arrow,  left)  and  two‐photon  (red  arrow,  right)  initiated transitions between molecular energy  levels. Z defines direction of  the  laser  beam.  Regions  that  are  appreciably  excited  are  shown  in  green.  2PA‐related processes are more localized in space.   Two‐photon  absorption  (2PA)  is  a  nonlinear  optical  process  in  which  chromophores,  that  would  normally  be  excited by a single photon (1PA) of ultraviolet light are excited  by  two  photons  of  near‐infrared wavelengths  (Figure  1).  The  probability  of  2PA  is  quadratically  dependent  on  intensity  of  the  incident  radiation,  which  opens  a  possibility  of  spatial  localization. A longer wavelength provides a better penetration  in  absorbing  or  scattering  media,  and  the  quadratic  dependence of the absorption probability on the light intensity  allows for high 3D spatial selectivity through the use of a tightly  focused  laser  beam.  Under  tight  focusing  conditions,  the  absorption  is  confined  at  the  focus  to  a  volume  of  order  3  (where  is the laser wavelength), and any subsequent process,  such  as  fluorescence  or  photoinduced  chemical  reactions,  is  localized in this small volume.  A  number  of  molecular  structures  were  designed    as  candidates  for  2PA  initiated  photochromic materials  [7],    but  were  found  to  have  insufficient  2PA  cross‐sections.  A  functioning 2PA based 3D optical storage system was reported  recently  [8].  It  offered  an  interlayer  separation  of  5  µm  and  depth  of  up  to  600  µm.  However,  the  diode  laser  commonly  used in DVD devices was not able to provide sufficient intensity  to record information, and a higher‐power solid‐state laser was  used  for  recording.  Thus,  design  of  photochromic  dyes  with  higher  2PA  cross‐sections  and  quantum  yields  is  vitally  important  for  practical  applications.  The  engineering  of  such  materials  is  in  need  of  rational  strategies  able  to  guide  the  progress in the material design.  There are  two aspects of  this  rational  strategy: prediction  of the two‐photon absorption cross‐sections of a candidate, and  prediction  of  its  photochemical  reactivity.  These  aspects  are  described  in  the  following  Subsections  2  and  3.  Each  of  them  includes  a  theoretical  background,  benchmark  study  and  relevant  applications.  The  final  Subsection  suggests  a  combination  of  methods  for  the  rational  design  of  a  2PA  initiated photochromic switch.  Computational  photochemistry  offers  a  number  of  theoretical methods for investigation of photochemical reaction  mechanisms.  Unlike  thermally  activated  chemical  reactions,  which  take  place  in  the  ground  electronic  state  (S0),  a  photochemical process involves the electronically excited state  (S1).  During  this  process  the  reactive  system  is  electronically  excited  from  S0  to  S1,  and  after  some  evolution  on  the  upper  potential energy surface (PES) decays back to the ground state  in  either  product  or  reactant  basin  through  conical  intersections (CIX) [9].  While accurate quantum dynamics calculations [10‐12] for  triatomic  systems  with  CIX  have  been  performed,  large  molecules  of  practical  interest,  necessitate  the  semi‐classical  [12,13]  approach.  In  this  approach  adiabatic  surfaces  are  treated  classically,  while  transitions  between  them  estimated  using  quantum  mechanics.  Another  approximation  involves  reducing  dimensionality  of  the  system  to  two  degrees  of  freedom  and  treating  the  resulting  two‐dimensional  problem  using  wave  packets  [14].  Further  simplification  is  called  Pathway  Approach  [15].  Instead  of  entire  PES,  it  considers  a  minimum energy path [16], which is followed by the center of  the  wavepacket  [17].  This  approach  is  focused  on  local  properties of  the PES, such as minima, barriers, and slopes.  In  this Pathway approach, CIX  serves as a  funnel, which delivers  the  excited  state  intermediate  to  the  ground  state  reactant  or  the product, so that quantum yield is largely determined there  [18]. It has been found very useful for qualitative analysis of the  reaction  mechanisms,  prediction  of  photoproducts,  and  rationalization of experimental excited state lifetimes, quantum  yields, absorption and emission spectra [19,20].  A  number  of  computational  tools  have  been  developed  to  predict  the  PES  [18,21,22].  Most  of  them  are  based  on  the  Complete Active Space (CAS) method, where molecular orbitals  are determined  self‐consistently within  a multiconfigurational  ansatz  (full  configuration  interactions  in  a  limited  orbital  space).  While  this  method  yields  acceptable  molecular  geometries,  the  account  for  dynamic  electron  correlation  is  required  for  accurate  energy  predictions,  and  is  typically  accomplished  by  the  multireference  configuration  interaction  or  perturbation  theory  treatment.  These  methods  are  computationally expensive and impractical for PES exploration  in large molecular systems.  Formalisms  based  on  density  functional  theory  (DFT)  are  being considered among the more affordable alternatives. One  of the approaches to treat excited states with DFT is the delta‐ Self  Consistent  Field  (ΔSCF), where  one  of  the  occupied  spin‐ orbitals is replaced with one of the unoccupied spin‐orbitals in  the  Kohn‐Sham  determinant,  and  orbital  relaxation  with  the  SCF  procedure  is  repeated.  In  general,  such  Kohn‐Sham  determinant  with  Fermi  hole  presents  a  saddle  point  and  special  orbital  optimization  algorithms  are  required  [23].  In  practice,  ΔSCF  is  possible  if  the  switched  orbitals  belong  to  different irreducible representations. This limitation is avoided  in  the  restricted  open‐shell  Kohn‐Sham  (ROKS)  method  [24],  which is based on the Ziegler sum rule [25], stating that a single  Kohn‐Sham  determinant  with  two  orbitals  occupied  by  unpaired electrons of  the opposite  spin  is  a 50/50 mixture of  pure singlet and triplet states and therefore has energy which  is  an  average  between  the  two.  This  method  was  shown  capable to describe geometry and perform molecular dynamics  of  the  open  shell  singlet  excited  state,  calculate  the  non‐ adiabatic coupling [26], and  include the S0–S1 surface hopping  into molecular dynamics simulations [27]. However, the double  excited state, describing the photoreactive system in pericyclic  minimum is not available in the ROKS methods.  An  alternative  approach  to  study  PES  of  electronically  excited systems is based on the time‐dependent (TD) or, more  precisely,  linear  response  DFT  formalism.  Instead  of  orbital  relaxation,  TD‐DFT  uses  an  alternative  procedure  where  the  Kohn‐Sham  (KS) wavefunction  is  expanded  in  terms  of  Slater  determinants, singly excited with respect to the reference state.  A rigorous  formulation of TD‐DFT [28] demonstrates  that  this  description  is  in  principle  exact,  given  that  the  frequency‐ dependent  exchange‐correlation  functional  is  known.  In most  practical  applications,  however,  this  frequency  dependence  is  ignored  (so  called  adiabatic  TD‐DFT).  This method was  often  reported  to  accurately  predict  electronic  spectra  and  excited  144  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161    state geometries. However, TD‐DFT was found to be somewhat  less  successful  in  description  of  PESs  in  the  vicinity  of  a  CIX  [29].  The  recent  study    [30]    demonstrated  that  these  difficulties may be routed in the failure of the restricted Kohn‐ Sham  formalism  for  the  reference  ground  state  close  to  the  geometry of the pericyclic minimum, and introduced a possible  solution, reviewed below in Section 3.    2. Prediction of two­photon absorption cross­sections  2.1. TD­DFT and a posteriori Tamm­Dancoff approximation  to the second­order TD­DFT  Density functional theory in the KS approximation [31] was  a method of choice  in  the solid state  theory  for a  long time.  It  received  recognition  as  a  reasonably  accurate  first  principles  method for the ground states of  large molecular systems after  the Generalized Gradient approximation was combined with a  fraction  of  Hartree‐Fock  (HF)  exchange  (hybrid  GGA)  [32].  Instead of a multiconfigurational wavefunction, the KS method  accounts  for  electron  correlation  through  the  exchange‐ correlation potential. In particular, contributions of the double  and  higher  excited  configurations  to  the  ground  states  are  included  in  the  single  Slater  determinant  KS  description  implicitly.  It  was  shown  that  the  exact  electron  density  obtained  in  multireference  ab  initio  methods  (such  as  Full  Configuration  Interaction)  can  be  mapped  onto  an  effective  single‐particle KS description  even  for molecules  far  from  the  equilibrium [33]. The ground state wavefunction is simply not  available  in  DFT  for  analysis,  and  the  single  KS  determinant  describes hypothetical system of noninteracting electrons used  in  the  KS  theory,  rather  than  a  molecular  system  of  interest.  Despite the single determinant appearance, DFT is not a single  reference method. This statement becomes apparent when the  Fermi  broadening  is  applied  to  the  occupation  numbers,  and  these  fractionally  occupied  KS  orbitals  are  compared  with  Natural  orbitals,  obtained  in  one  of  wavefunction  based  methods [34]. From a practical standpoint, approximate hybrid  exchange‐correlation functionals (such as B3LYP [35]) result in  energies close to coupled cluster models in quality [32]. There  are  several  extensions  of  DFT,  developed  to  describe  the  excited states. One approach is to use multireference strategies  [36‐39]. More often, however, DFT  is  combined with  the  time  dependent  perturbation  theory  approach,  which  we  shall  briefly review next.  Instead  of  obtaining  lowest  energy  solutions  to  the  Schrödinger  equation  based  on  the  variational  principle,  excited  states  can  be  described  by  the  time‐dependent  perturbation  theory  treatment  as  interaction  between  the  molecular system and the electric field [40]. Several formalisms  had been developed to obtain excitation energies and transition  dipoles  directly  [41],  including  the  equation  of motion  (EOM)  and  polarization  propagator  approaches.  When  the  terms  of  second  and  higher  order  in  the  external  filed  are  neglected  (linear  response,  LR),  and  the  HF  method  is  adopted  for  description  of  the  ground  state,  the  approach  is  called  Time‐ Dependent  Hartree‐Fock  (TDHF)  method,  also  known  as  random  phase  approximation  (RPA).  When  applied  to  DFT  ground state [28,42],  the LR approximation  is often called LR‐ DFT, or TD‐DFT. The LR‐DFT method results in non‐Hermitian  eigenvalue problem:                       Y X Y X AB BA                                                             (1)    solution to which are the excitation energies Ωα and transition  density matrices  ξα  for  the ground  to excited  state  transitions.  In the basis of occupied (i,j) and vacant (a,b) KS orbitals of σ,τ   subsets  (σ,τ=α,β),  transition  density  is  block‐diagonal  with  occupied‐vacant  iaX )(   and  vacant‐occupied  aiY )(   blocks being nonzero. Matrices A and B are defined as      bjaiiaijabbjai KA ,, )(      and      jbaibjai KB ,,                                                                             (2)    For the hybrid DFT with cHF fraction of the HF exchange, the  coupling matrix K is expressed through the second derivatives  of  the  exchange‐correlation  functional  w,  Coulomb  and  exchange integrals as:    )|()|()||)(1(, ijabcjbiajbwiacK HFHFbjai                   (3)    The  RPA,  or  TDHF,  equation  is  then  a  limiting  case  with  cHF=1, where the matrix A consists of interactions between two  singly‐excited  configurations  )||( jbHia  ,  also  known  as  the  configuration  interaction  limited  to  singles  (CIS)  Hamiltonian.  The  matrix  B  includes,  by  virtue  of  swapping  indexes,  the  excitations  from  virtual  to  occupied  molecular  orbitals  (deexcitations)  of  the  form  )||( bjHia  .  Mathematically,  they  are  equivalent  to  the  matrix  elements  between  the ground and doubly excited states  [43]. Thus, LR‐ DFT  accounts  for  double  excitations  in  two  ways:  implicitly,  though  the  exchange‐correlation  functional,  and  explicitly,  through the deexcitation matrix B.  Sometimes  an  additional  approximation  is  introduced  in  the  LR  formalism.  It  is  called  Tamm‐Dancoff  approximation  (TDA)  [44,45],  and  consists  in  neglecting  the  deexcitation  matrix B in equation (1). When applied to the HF ground state,  it  is equivalent to the CIS method. Double excitation character  is  included  in  the TDA‐DFT  formalism only  implicitly  through  the approximate exchange‐correlation (XC) potential. TDA was  found to be accurate when an exact XC potential was restored  from an exact electron density in systems close to equilibrium  [46]  (but  not  when  covalent  bond  is  stretched  [47]).  To  a  certain  extent,  TDA  may  correct  the  deficiencies  in  the  approximate exchange‐correlation functionals [48].  In practice,  an  implicit  account  for double  excitations was  found  to  be more  important  than  the  role  of  the  deexcitation  matrix. For instance, the study of excited states in radicals [49],  demonstrated that TD‐DFT gives results comparable to TD‐HF  for single excited states and more accurate excitation energies  for  the  double  excited  states.  The  results  for  linear  polyenes  were mixed. Several studies  [48,50,51] had reported a correct  state  ordering  only  with  pure  XC  functionals  (and  not  with  more  accurate,  otherwise  hybrid,  functionals),  combined with  Tamm‐Dancoff approximation and diffuse basis functions. Both  successes  and  failures  of  LR‐DFT  to describe potential  energy  surfaces  of  single  and  double  excited  states,  conical  intersections  between  states,  and  photochemical  transformations have been reported [29,52,53].  The  alternative  approaches,  taking  one  lowest  double‐ excited  states  into  account  explicitly  at  the  LR  level  had  been  recently  developed.  Spin‐flip  DFT  [54],  and  noncollinear  XC  with spin‐flip excitations [55] use the (HOMO)1(LUMO)1 triplet  (where  HOMO  and  LUMO  stand  for  the  highest  occupied  and  lowest  unoccupied  molecular  orbitals  respectively)  as  a  reference state and obtain both the ground and double excited  states  as  single  excitations.  Another  approach,  called  Dressed  TD‐DFT  [56]  introduces  frequency  dependence  into  the  XC  kernel  (non‐adiabatic  approximation  to  TD‐DFT)  by means  of  adding  matrix  elements  involving  double  excited  HOMO2  →LUMO2 configurations  into matrices A and B of eq.  (1). As a  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  145  result, the 2Ag state acquires the double excited character, and  its energy is considerably lowered. Although double excitations  could be included in extended or higher‐order RPA formalisms  [57], these methods did not receive wide attention. Instead, LR  has been applied to correlated ground states,  including  Multi‐ Configuration Self Consistent Field (MCSCF) [58], and Coupled  Cluster (CC)  [59‐61] methods.  Let us focus on the transition dipole moments next.  In LR‐ DFT  they  are  readily  obtained  as  convolution  of  the  dipole  moment operator with the transition densities:     )(   Tr                                                                              (4)    State‐to‐state transition dipoles μα,β do not appear in the LR  approximation,  unless  the  excited  state  α  is  taken  as  the  reference state. In order to obtain the expressions for μα,β one  has to extend time dependent perturbation theory to Quadratic  Response  (QR). This extension was  initially developed  [62,63]  for HF and MCSCF reference states using an explicit exponential  unitary  time  dependent  transformation.  The  results  of  the  perturbation  treatment  were  expressed  in  terms  of  response  functions  (first,  second,  and  third‐order  corrections  to  the  expectation  values  of  an  arbitrary  operator). Many molecular  properties may be extracted from single and double residues of  these  response  functions  [64]  at  the  resonant  frequencies  (poles of the response function). Specifically, a transition dipole  between  excited  states  can  be  expressed  through  the  ground  state  dipole  moment  μo,o,  ground‐to  excited  state  dipole  moments  μα  and  μβ,  and  the  second  residue  of  the  dipole  quadratic  response  function  with  two  electric  dipole  perturbations μb and μc  [63]:    )](,;)(lim[lim)( 1 0,0,       b cba c cbaa cb      (5)    Substitution of  the LR values  in place of dipoles evaluated  with the exact states leads to summation over large number of  states.  The  explicit  summation  can  be  replaced  by  iterative  solution of the linear equations [63], which may be recast [65]  in a  form, similar  to eq.  (1). Luo et al.  studied  linear polyenes  (butadiene,  hexatriene,  and  octatetraene)  at  the  QR‐HF  level  [66].  They  found  that  the  1Bu  state  dominates  the  linear  spectra, while  three or  four Ag  states  (depending on  the basis  set),  that  have  large  transition  dipoles  from  the  1Bu  state,  appear  on  the  two‐photon  absorption  spectra  with  maxima  near 1.5–1.7 of the bandgap (excitation energy of the 1Bu state).  The  quadratic  response  theory  combined  with  DFT  was  also  used  for  simulations  of  two‐photon  absorption  spectra  [67].  The  quadratic  response  formalism,  extended  to  CC  [68]  and  DFT  [69‐71]  reference  states  is  implemented  in  the  Dalton  program [72].   An  alternative  formulation  of  time‐dependent  perturbation  theory  for  excited  states  is  known  as  Coupled  Electronic  Oscillator  (CEO)  approach  [73,74].  It  uses  density  matrix  (Liouville  space)  representation  and  is  based  on  the  classical  Heisenberg  equation  of  motion  for  the  ground  state  density. It was recently extended from the HF to DFT reference  state  [75]. When  only  terms  of  the  first  order  in  the  external  field are retained, equation (1), equivalent to the LR formalism  is  obtained.  In  the  second  order  approximation,  the  solutions  are sought in the basis of LR solutions (transition densities). As  a  result,  linear  excitations  remain unchanged  in  the quadratic  formalism, and combined states ξβξα of a double excited nature  are added to the picture. The excitation energy for each of these  new states is equal to the sum of single excitations:      .                                                                          (6)    The second order CEO gives the transition dipole between  the ground and this double excited state as                    0 ,0 ))2((          VV ITr perm  (7)    Here  the  first  summation  runs  over  symmetrized  permutations  of  the  indexes,  second  summation  includes  negative  indexes,  I  is  the  identity matrix, ρ  is  the ground state  density matrix,  and Vαβ­γ  is  the  exchange‐correlation  coupling  term, expressed via KS operators V(ξ) on transition densities:      perm VITrV    ))()2(( 2 1                                      (8)    Further,  the  transition  dipole  between  the  double  excited  state and any other excited state is zero unless the other state  presents one of the components of this double excited state:      , ;     0,                                                                   (9)    The transition dipole between two single excited states is                     0, , ))2((          VV ITr perm  (10)    Thus, in the second order CEO the first double excited state  of  the  (HOMO)2→(LUMO)2  type  is  always  twice  higher  in  energy  than  the  HOMO→LUMO  excited  state  and  (unlike  in  wavefunction  based  methods),  it  never  mixes  with  single  excitations.  The  perturbative  treatment  can  be  further  extended to the third order, which results in corrections to the  states  and  leads  to  the  mixing  of  single,  double  and  triple  excitations within TD‐DFT formalism (eq. G5‐G7 in Ref. [75]):    , !3 1 0 0 ††† 0 )0(          gaaa V g                         (11)    , !2 1 0 0 †† 0 †)1(           gaa V ga                      (12)      , !2 1 2 !2 1 0 0 †††† 0 0 † 0 ††)2(                                 gaaa V a V ga V gaa        (13)    where  0 g ,  0 † ga ,  0 †† gaa  , and  0 ††† gaaa   stand for  the  ground,  single,  double,  and  triple  excited  states  of  the  uncoupled  system,  respectively.  An  alternative  perturbative  corrections  to  the  LR‐DFT  excitation  energies  were  proposed  recently by Ziegler et al. [76,77].  Double excited states do not appear when TDA is  invoked,  and  exchange‐correlation  coupling  terms  Vαβ­γ  vanish.  As  a  result,  state‐to‐state  transition  dipoles  coincide with  the  ones  obtained in the CIS method, defined by the first term of eq. (10).  Recently  we  suggested  [100]  to  apply  Tamm‐Dancoff  approximation  and  annihilate  Y  component  of  transition  density a posteriori,  after  LR  equation  had  been  solved.  Thus,  excitation  energies  and  ground  to  excited  transition  dipoles  remain unaffected, while CIS  formulas are applied  to calculate  146  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161    state‐to‐state  transition  dipoles  and  differences     between  the permanent excited and ground state dipole moments:     ))2((,    ITr ; ))2((    ITr ;         0 X        (14)  In addition, the double excited states    are introduced  with the following excitation energies and transition dipoles:      ; ))2((,0    ITr ;      , ;  0,                                                                     (15)    We  called  this  scheme a posteriori  Tamm‐Dancoff  (ATDA)  approximation  [100];  it  is  intermediate  between TDA and  full  second  order TD‐DFT. ATDA  allows  to  calculate  second  order  properties without solving equations of the full QR‐DFT, using a  simple modification to existing LR codes. Unlike QR‐DFT, ATDA  inherits  double  excited  states  from  second  order  CEO  formalism. In the Subsections 2.2‐2.3 we demonstrate accuracy  of  ATDA  in  calculation  of  transition  dipoles  between  excited  states  in  trans‐butadiene,  as  well  as  permanent  dipole  moments  for  the  excited  states  in  para‐nitroaniline,  and  compare the numerical values with the ones obtained at higher  theory  levels.  Comparison  of  the  one‐photon  and  two‐photon  absorption spectra  for some chromophores with experimental  data, made  in  the  Subsection  2.4,  illustrates  the  ability  of  the  ATDA to reproduce optical properties with a good accuracy.    2.2.  Benchmarking  ATDA  study  of  the  transition  dipole  moments   Conjugated  hydrocarbons  (also  known  as  polyenes)  and  their derivatives present an important class of compounds with  rich  photophysical  and  photochemical  properties.  These  properties  originate  in  a  highly  polarizable  π‐electron  system  and  find a wide use  in organic  electronics  [78]  and photonics  [79] applications. Another reason  for  interest  in electronically  excited states of polyenes is their role in biological processes of  vision and photosynthesis. Theoretical description of electronic  structure and electronic excited states in conjugated molecules  plays  a  critical  role  in  understanding  the  natural  and  engineered processes, and may assist  in  the rational design of  the new materials with improved properties.  Polyenes  often  served  as  a  testing  ground  for  new  theoretical  methods,  and  comprehensive  review  of  the  published results seems to be impossible. Presently consensus  is reached about ordering of the lowest excited states, with 2Ag  one being above 1Bu  (but close  in energy)  for  trans butadiene  and all‐trans hexatriene, and below the 1Bu state for the higher  hydrocarbons.  Transition  dipole  moments  between  excited  states of polyenes were much less studied, and double‐excited  nature  of  2Ag  states  is  still  a matter  of  discussion  [51]. At  the  same  time,  double‐excited  nature  of  selected  excited  states  remains  to  be  an  important  challenge  of  TD‐DFT  [47].  These  aspects of electronic structure, as well as the ability of TD‐DFT  methods  to  describe  them  present  one  of  the  focuses  of  this  review.  Accurate  numerical  values  of  transition  dipoles  between  excited  states  of  molecules  are  important  for  prediction  of  nonlinear  photonic  processes  [80,81],  such  as  excited  state  absorption  and  2PA.  States,  involved  in  these  processes,  are  often not observable  in  linear absorption spectra and some of  them  appear  to  possess  strong  double  excited  character  [82‐ 84].  The  relation  between  two‐photon  absorption  and  other  nonlinear optical properties of the ground sates is provided by  the  Sum  over  States  (SOS)  expression,  derived  within  the  perturbation  theory  approach  [85].  Several  authors  [86‐88]  had  noticed  that  first  and  second  hyperpolarizability  values,  obtained for the linear polyenes at a semiempirical theory level  within the SOS approach, are dominated by contributions from  a few, so called “essential” states, so that all other states can be  excluded  from  consideration.  This  approximation  appeared  attractive for structure‐property relationship schemes [89,90],  but  later  ab  initio  studies  [66,91]  reported  it  to  be  an  oversimplification. In this Subsection we review the results of a  more  advanced  correlation  treatment  [92],  and  identify  approximations  responsible  for  disagreements  between  different ab initio approaches.  In this Section we aim to  focus on the valence ππ*‐excited  states,  which  are  primarily  responsible  for  the  second  and  higher  polarizabilities  of  conjugated  hydrocarbons.  Linear  all‐ trans  polyenes are planar molecules of C2h  symmetry,  their π‐ orbitals  belong  to  au  and  bg  irreducible  representations.  Configurations  (Slater  determinants)  describing  the  electron  transitions  between  the  orbitals  of  the  same  symmetry  contribute to Ag states, while transitions between the orbitals of  different  symmetry  contribute  to  Bu  states.  Since  the  ground  state  is  1Ag,  one‐photon  transitions  to  Bu  states  are  dipole‐ allowed and to Ag states are dipole‐forbidden in the one‐photon  regime, while two‐photon transitions are forbidden to Bu states  and allowed to Ag states. Besides the spatial symmetry, π‐states  can  be  classified  according  to  the  so‐called  alternacy,  or  particle‐hole symmetry  [93].  It  is exact only with some model  Hamiltonians  (Huckel  or  Parriser‐Parr‐Pople),  and  becomes  approximate  after  σ‐π  and  second‐neighbor  interactions  are  included. However,  it  is  useful  for  interpretation  of  transition  dipoles  from  ab  initio  calculations.  Let  us  first  introduce  shorthand notations 1, 2, …, m for the occupied and 1’, 2’, …, m’  for the vacant π‐orbitals, with HOMO being 1, and LUMO being  1’.  When  model  Hamiltonians  are  used,  configurations m→n’  and n→m’  are degenerate  and  therefore  contribute  equally  to  the  state  wavefunction.  Their  in‐phase  and  anti‐phase  linear  combinations  give  rise  to  plus  and  minus  states  [94]  The  valence  bond  theory  description  of  in‐phase  states  [95]  includes  a  substantial  contribution  from  ionic  resonance  structures, while  for  the anti‐phase states  it does not. For  this  reason plus and minus states are often called ionic and covalent  respectively  [96].  The  actual  sign  of  the  amplitude  for  every  given  configuration  depends  on  the  signs  of  singly  occupied  orbitals included in this configuration and is assigned arbitrary  in  most  of  software  packages.  Hence,  we  will  use  ionic  and  covalent  notations,  as  plus  and  minus  notations  tend  to  be  misleading.  Similar  to  the  spatial  symmetry,  the  alternacy  symmetry results in selection rules for one‐photon transitions,  as  the dipole  transition moment between any  two  ionic  states  or  between  any  two  covalent  states  is  zero.  The  ground  state  1Ag behaves like a covalent state. The excited configurations of  the type of m→m‘ behave like ionic states for singlet‐spin states  and  covalent  states  for  triplet‐spin  states.  The  doubly  excited  configurations of mm→n’n’ type behave like covalent states and  mix with the singly excited covalent configurations. This results  in multiconfigurational character of  the covalent 2Ag  state and  explains why this state is strongly stabilized, sometimes below  the ionic 1Bu state, of predominantly 1→1’ nature.  Ab­initio  description  of  the  molecular  excited  states  typically  starts  with  molecular  orbitals,  optimized  in  the  self  consistent  field  (SCF)  procedure  with  a  single  Slater  determinant  wavefunction  of  the  ground  state,  known  as  the  restricted  HF  method.  Excited  states  can  in  principle  be  described by single determinants (or a couple of determinants,  for open shell singlets), which are constrained to be orthogonal  to  Slater  determinants  of  the  lower‐lying  states.  While  this  approach works  reasonably well  at  the DFT  level  [26,97],  the  HF description of excited states is generally inaccurate. Instead,  configuration interaction (CI) methods are used, where one or  more  orbitals  in  the  HF  determinant  are  substituted  with  unoccupied  orbitals  to  form  excited  configurations.  The  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  147  wavefunction  of  the  system  is  expressed  as  a  linear  combination  of  these  configurations  and  electronic  states  are  found  by  diagonalization  of  the  Hamiltonian.  All  possible  substitutions in the HF determinant yield the Full Configuration  Interactions  method  (FCI).  It  gives  an  exact  solution  to  the  Schrödinger  equation  for  a  given  atomic  basis  [98].  The  exponential growth of the computational effort with the size of  the system makes FCI attainable only for very small molecular  systems, and  for practical reasons various  truncation schemes  are introduced. The simplest scheme restricts the expansion of  the  wavefunction  to  single  substitutions  (CI  Singles,  or  CIS),  another  one  limits  the  substitutions  to  singles  and  doubles  (CISD).  To  reduce  the  computational  effort,  amplitudes  for  double  substitutions  can  be  evaluated  perturbatively  (instead  of  the  variational  approach),  resulting  in  the  CIS(D)  method  [99]. When applied to the ground state, the method is known as  MP2  (Moeller‐Plesset  second  order  perturbation  theory).   Second‐order  Algebraic  Diagrammatic  Construction  presents  yet  another  perturbation  correction  method,  where  the  correction is applied to the matrix elements before solving CIS  equations  [100].  Variational  CISD  description  for  the  ground  state  is  known  not  to  be  size‐consistent,  which  decreases  its  accuracy  for  the  larger  molecules.  In  order  to  maintain  size‐ consistency,  selected  triple  and  quadruple  excitations  are  added to the wavefunction in the form of products of single and  double excited configurations. Their amplitudes are also taken  as  products  rather  than  evaluated  independently.  This  approach  results  in  coupled  cluster  expansions  (CCSD,  CC2,  etc.)  [101].  At  the  moment  couple  cluster  methodology  provides  arguably  the  best  accuracy/cost  ratio  among  wavefunction based methods. The CI method, based on the CC  reference  (SAC‐CI)  [102,103]  extends  this  treatment  to  the  excited  states. Both  single  and double  (SD‐R)  [104]  and up  to  sixth‐order  (General‐R)  [103]  excitation  truncation  schemes  were developed. Both SD‐R and General‐R were shown to yield  similar results for the 2Ag state of trans‐butadiene.  The  results  for  the  singlet  ππ*‐states  of  butadiene  are  presented  in  Table  1.  The  reader  can  find  these  and  other  results  on  hexatriene  and  octatetraene  in  Ref.  [92].  The  chemical structure of trans butadiene is shown in Scheme 2.        Scheme 2. Trans‐butadiene.    A close examination of the major configurations in the SAC‐ CI  wavefunction  reveals  that  the  alternacy  symmetry  classification  into  ionic  and  covalent  states  approximately  holds.  Out  of  two  linear  combinations  of  (i­j’)  and  (j­i’)  determinants, the covalent combination is lower in energy and  contains a larger contribution from the double excitations. The  actual  signs  of  the  amplitudes  depend  on  the  phase  of  molecular orbitals, and are therefore arbitrary. For covalent Ag  states, one can observe a large transition dipole from the ionic  1Bu  state.  For  covalent  Bu  states,  transition  dipoles  from  the  covalent ground state are much lower,  than for  ionic Bu states,  according to alternacy selection rules. As one can see from SAC‐ CI  wavefunction  analysis,  this  stabilization  of  the  covalent  vs.  ionic state corresponds to 7‐14% contribution from the leading  double  excited  configuration,  which  is  found  insignificant  in  ionic states. The same stabilization is paralleled by both CIS(D)  and TD‐DFT excitation energies. For instance, in hexatriene the  CIS  covalent  2Bu  state  is  0.5eV  above  the  ionic  3Bu  state,  the  CIS(D)  covalent  2Bu  state  is  1.5  eV  below  the  ionic  3Bu  state.  Similarly,  the  ATDA  covalent  2Bu  state  is  1  eV  lower  than  the  ionic 3Bu state. While this stabilization is achieved through the  perturbative correction for double excitations in case of CIS(D),  it is afforded by means of the exchange‐correlation potential in  the case of TD‐DFT. Comparison SAC‐CI with CIS results, where  covalent  states  are  systematically  destabilized with  respect  to  ionic states, is especially instructive.   Thus, analysis of SAC‐CI wavefunction, which includes both  static  and  dynamic  electron  correlations  shows  that  double‐ excited character of 2Ag and other lower‐lying covalent states in  the linear polyenes is relatively weak (less than 20%) and these  states  can  be  therefore  classified  as  predominantly  single  excitations. This is in contrast with ionic states that have purely  single  excited  character.  Some  of  the  higher‐lying  states  however have considerably higher contributions of the doubly  excited  configurations  (40‐80%).  Excitation  energies  of  these  states  are  close  to  the  sum  of  the  excitation  energies  of  the  respective  predominantly  single  excited  states.  Even  though  doubly  excited  configurations  appear  to  contribute much  less  than 100% to the wavefunction, this contribution is more than  twice  greater  than  the  one  observed  in  predominantly  single  excited  states.  Thus,  these  states  can  be  interpreted  as  predominantly doubly excited. This interpretation is confirmed  by  the  trends  in  excitation  energy  and  transition  dipoles,  as  described below. Thus,  single  and double  excitations  in  linear  polyenes at the equilibrium geometry are separated by a gap in  double  excited  character  (below  20%  and  above  40%).  It  is  conceivable that far from equilibrium this observation may not  hold.  Let  us  consider  the  4Ag  state  of  butadiene  as  an  example.  According to Table 1, 4Ag has excitation energy of 12.6 eV, close  to  the  double  excitation  energy  for  the  1Bu  state  (6.8  eV).  Meanwhile,  its  leading  configuration  (11‐1'1')  corresponds  to  the double HOMO‐LUMO excitation (1‐1'), which is the leading  configuration  in  the  1Bu  state.  Similarly,  the  3Bu  state  of  butadiene with the leading configurations of (11‐1'2')‐(12‐1'1')  originates from the product of singly excited 1Bu (1‐1') and 3Ag  (1‐2')‐(2‐1')  states.  Its  excitation  energy  of  15.3  eV  is  nearly  equal to the sum of the respective excitation energies of 6.8 and  9.4 eV.  The  relationships  observed  in  SAC‐CI  results  are  closely  resembled  by  the  ATDA‐B3LYP  level  of  theory.  The  2Ag  and  other  states,  appeared  to  be  predominantly  single  excitations  with SAC‐CI, were found to have the excitation energies closer  to  benchmark  values  than  uncorrected  SAC‐CI/DZp  values.  Unlike CIS(D) and SAC‐CI results, covalent states are stabilized  relative  to  the  corresponding  ionic  states  without  explicitly  including  double  excited  configurations  into  the  description.  This confirms that the same electron correlation effects, which  are  described  by  double  excited  configurations  in  the  wavefunction, are accounted for by the approximate exchange‐ correlation  functional  B3LYP.  On  the  other  hand  the  states,  interpreted  as  predominantly  doubly  excited,  are  well  reproduced  by  the  second  order  CEO  extension  to  TD‐DFT,  after  eqs.  (6–10)  are  simplified  with  ATDA.  These  states  are  completely  absent  in  the  list  of  CIS  and  CIS(D)  states.  It  is  necessary  to stress,  that double excited states are not present  in  the  output  of  standard  computer  codes,  implementing  TD‐ DFT in either linear or quadratic response.  Transition dipole moments  calculated with ATDA‐DFT are  in  better  agreement  with  SAC‐CI  values  than  CIS  ones.  They  also  compare  favorably with QR‐DFT results. While ground  to  Bu  transition  dipole moments  are  the  same  at  ATDA  and  QR‐ DFT methods  (being  LR  values),  1Bu  to  nAg  transition  dipoles  are different. Figure 2 illustrates the comparison of 1Bu  to nAg  transition dipoles in butadiene obtained in the ATDA, QR‐DFT,  and  SAC‐CI.  Ionic  states  are  almost  equally  dark with QR  and  ATDA. Covalent Ag states of mostly single excited character are  bright  in SAC‐CI, darker with ATDA, and even darker with QR.  Covalent  states  of  mostly  double  excited  character  are  completely  absent  in  QR,  but  their  transition  dipoles  are  borrowed  by  other  states  close  in  energy.  Even  with  this  in  mind, ATDA values are still closer to SAC‐CI ones.   148  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161      Table 1. Excitation energies ∆E (eV), major contributing configurations, and % contributions of the leading double excited configurations for the lowest π‐states  of trans butadiene. Transition dipoles μ (a.u.), from the ground state to nBu states, and from 1Bu to nAg states are also reported. All basis sets are DZp. The  benchmark values of excitation energies obtained at CAS‐PT2/ANO are taken from Ref. [105].      State  1Bu  2Bu  3Bu  2Ag  3Ag  4Ag  5Ag  CIS  ∆E,eV  6.66  12.15  ‐ 9.66 9.03 ‐  ‐ Config.  (1‐1')  (2‐2')  ‐ (1‐2')+(2‐1') (1‐2')‐(2‐1') ‐  ‐ μ,a.u.  2.661  0.784  ‐ 2.682 0.112 ‐  ‐ CIS(D)  ∆E,eV  6.77  11.33  ‐ 7.95 9.28 ‐  ‐ SAC­CI (SD­R)  ∆E,eV  6.86  11.72  15.31 7.51 9.41 12.65  16.71 Config.  (1‐1')  (2‐2')+(11‐1’2')  (11‐1'2’)‐(12‐1’1')  (2‐1')+(1‐2')+(11‐1’1')  (1‐2')‐(2‐1') (2‐1')+(1‐2')‐  (11‐1’1')  (12‐1’2’)+(22‐1’1’)+  (11‐2’2’)  %D  0  2  78 13 0 53  83 μ,a.u.  2.291  0.874  0.17 1.695 0.306 2.698  0.557 ATDA­B3LYP  ∆E,eV  5.94  10.33  14.52 7.11 8.57 11.89  16.28 Config.  (1‐1')  (2‐2')  (1‐1')((1‐2’)‐(2‐1’)) (2‐1')+(1‐2') (1‐2')‐(2‐1') (1‐1')(1‐1')  (1‐1’)(2‐2’) μ,a.u.  2.131  0.706  1.485 0.086 2.131  0.706 QR­B3LYP  μ,a.u.      1.293 0.199 ‐  ‐ CAS­PT2/ /ANO  ∆E,eV  6.06  ‐  ‐ 6.27 ‐ ‐  ‐   According  to  eq.  (11–13),  transition  dipole  moment  between  the  single  and  the  double  excited  state  (including  this  singly  excited state as its component) is equal to one from the ground  to  the  second  component  of  the  double  excited  state.  SAC‐CI  transition dipole values between respective states are close to  this  CEO  relationship,  which  confirms  our  classification  of  covalent  states  into  predominantly  single  and  predominantly  double  excitations.  The  accuracy of ATDA  in prediction of  the  permanent dipole moments is reviewed in the next Subsection.      Figure 2. Calculated transition dipoles from 1Bu state to a manifold of nAg states  in  trans‐butadiene  vs.  excitation  energy  of  nAg  states.  With  text  labeled are configurations and percentage of double excitation. Color of the  text  correspond  to  the  legend:  red  is  for  SAC‐CI,  and  blue  is  for  ATDA.  Circled are the states of similar configuration in different methods. The first  two  Ag  states  below  10eV,  obtained  in  the  ATDA,  are  of  single  excited  nature. All other ATDA states above 10eV with essentially non‐zero dipoles  correspond to double exited states constructed as (1Bu nBu).    2.3.  Benchmarking  ATDA  study  of  the  permanent  excited  state dipoles  In  order  to  investigate  the  permanent  dipole moments  of  the  excited  states  in  a  prototypical  donor/acceptor  (D/A)  substituted  conjugated  molecule,  let  us  consider  para‐ nitroaniline (PNA). Its structure is shown in Scheme 2.      Scheme 2. Para‐nitroaniline molecule.    PNA has a nonbonded electron pair localized on the amino  group  (which  plays  a  role  of  donor)  and  a  low‐lying  vacant  antibonding  orbital  localized  on  the  nitro  group  (electron  acceptor). The ‐conjugation through the six‐membered phenyl  ring allows for partial transfer of the electron density from the  donor  to  the  acceptor  groups  in  both  the  ground  and  excited  states. This electron mobility results in a number of properties  important  for  organic  electronics  [78]  and  photonics  [79]  applications.  Photoinduced  change  in  the  permanent  electric  dipole  moment  reflects  the  degree  of  charge  transfer  resulting  from  electronic excitation [106], and is necessary for understanding  of many practically important systems, from isolated molecules  in  the  gas  phase  and  in  solution,  to  the  band  structures  of  organic  semiconductors  and  crystalline materials,  to  complex  biological systems [107]. Knowledge of the excited state dipole  moments of PNA homologs  is also    important  for wide variety  of light induced charge transfer reactions [108].  Accurate  numerical  values  of  permanent  dipole  moments  (as  well  as  transition  dipole  moments)  are  also  useful,  along  with  transition  dipole  moments,  for  prediction  of  nonlinear  optical  properties  with  the  SOS  formalism  [85].  The  reliable  estimation  of  nonlinear  optical  properties  typically  requires  higher  theory  levels  that  include  electron  correlation.  Prediction  the  dipole  moments  of  the  excited  states  was  a  challenge  for  TD‐DFT  [47].  However,  in  the  previous  Section,  we have shown that ATDA to the second‐order CEO formalism,  applied  to  the  TD‐DFT  can  be  used  for  an  accurate,  yet  computationally  inexpensive  prediction  of  the  transition  dipoles. The focus of this Section is to analyze a performance of  the ATDA method in prediction of permanent dipole moments  of excited states based on comparison with the results obtained  at a high‐level wavefunction theory. As a test case, we analyze  electronic structure of the five lowest excited states in each of  four  irreducible  representation  of  the  C2h  symmetry  group  of  para‐nitroaniline.  We  optimized  the  geometry  of  para‐nitroaniline  in  assumption of C2v symmetry, at  the CCSD/DZp  level of  theory.  The B3LYP exchange‐correlation functional was used in all DFT  calculations.  The  basis  set  DZp,  due  to  Nakayama  et al.  [109]  was  obtained  by  deleting  p‐functions  on  hydrogen  from  the  correlation‐consistent cc‐pVDZ basis set [110] and used for all  the  calculations  in  this  Section.  The  properties  of  the  excited  states in the EOM‐CCSD method were obtained with the ACES2  program [111]. The permanent dipole moments at the QR‐DFT  level  were  obtained  with  the  Dalton  2.0  [72].  The  Gaussian  2003  suite  of  programs  (revision  E.1)  [112]  was  used  to  perform  TD‐DFT  calculations  in  two  different  ways.  The  numerical  derivatives,  according  to  Eq.  (1)  were  calculated  using the keyword ‘Polar’. The calculations with the unrelaxed  TD‐DFT  density  were  also  carried  out  using  the  keyword  ‘Density=RhoCI’.  A  locally modified  version  of  the  Gaussian98  [113] was  used  for  the  calculations  at  the  ATDA‐DFT  level  of  theory.  Properties of the five lowest singlet excited states for each  irreducible  representation  within  the  C2v  symmetry  are  collected  in  Table  2.  The  properties  include  the  excitation  energies  (E),  the  projection  of  permanent  dipole  moments,  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  149  which  coincides with  the  rotation  axis  in  the  C2v  point  group,  and  leading  electronic  configurations  of  the  excited  states  obtained  at  EOM‐CCSD  and  TD‐DFT  theory  levels.  The  EOM‐ CCSD  excitations  in  each  symmetry  species  are  sorted  by  energies.  For  TD‐DFT  method  the  states  are  ordered  by  the  close  correspondence  of  their  leading  configurations  to  the  leading configurations of the EOM‐CCSD excited states. In order  to  describe  configurations,  we  use  the  same  shorthand  notations  for  the  occupied  and  vacant  orbitals  in  each  symmetry  species  as  in  the  previous  Section.  The  irreducible  representations of the C2v symmetry are assigned the following  numbers: (1) a1, (2) b2, (3) b1, (4) a2. The number of irreducible  representation  is  specified  in  brackets  after  the  orbital.  For  instance, HOMO in the b1 symmetry is labeled 1(3), LUMO in the  b2 symmetry is labeled 1’(2), and a configuration obtained from  the HF determinant by excitation of an electron from one  into  the other is denoted (1(3)‐1'(2)). Close  examination  of  the  EOM  excited  states  reveals  that  contribution of double‐excited determinants to the lower‐lying  states  of  para‐nitroaniline  is  relatively  small  (less  than  4%).  This  is  typical  for  the  chromophores  of  the  donor‐acceptor  type. All four methods of calculating the dipole moments within  the  TD‐DFT  formalism  report  different  values  for  the  same  states. While Quadratic Response values are reasonably close to  the  numerical  differentiation  results,  the  permanent  dipole  moments  for  the  TD‐DFT  obtained  with  the  unrelaxed  electronic  density  (‘unrelaxed  dipoles’)  are  often  overestimated. ATDA dipoles are in better agreement with the  EOM  results,  especially  for  the  lowest  excitations  in  each  irreducible representation. The  deviations  of  QR‐DFT  dipoles  from  the  numerical  derivatives  for  some  states  may  be  explained  by  a  large  discrepancy  between  QR‐DFT  and  EOM‐CCSD  energies.  The  second  order  corrections  for  excitation  energies  are  not  calculated  in  the  QR‐DFT  implemented  in  the  Dalton  2.0.  Permanent dipole moments are obtained  in  the QR  formalism  as  a  projection  of  the  response  function  onto  the  energy  spectrum calculated in the LR. For those states whose energies  are  not well‐defined  in  the  LR,  i.e.  their  second  order  energy  corrections  are  going  to  be  large,  the  discrepancies  in  dipole  moment values can be big. That deviation can be observed for  near‐degenerated levels, so that if energies of two states differ  less  that  0.1eV,  they  interfere  in  the  QR  formalism.  Another  possible explanation is switching the state ordering in the finite  field  used  for  numerical  differentiation  that  can  lead  to  incorrect calculation of the derivatives.  In  order  to  quantify  the  differences  between  EOM‐CCSD  and  TD‐DFT  excitation  energies  and  permanent  diploes,  we  calculated  the  standard  deviations.  They  are  reported  in  the  last two rows of Table 2  One can see that numerical derivatives, QR‐DFT and ATDA‐ DFT  values  are  close  to  each  other  and  to  coupled  cluster  values,  while  unrelaxed  dipole moments  deviations  are  three  times  as  large.  We  find  ATDA‐DFT,  QR‐DFT  and  numerical  energy  derivatives  to  be  in  good  agreement  with  the  results  obtained  with  the  EOM‐CC  approach.  Unrelaxed  density  TD‐ DFT calculations, on the other hand, yield numerically  inferior  results. This validates the accuracy of ATDA for the description  of  permanent  dipole  moments  of  excited  states  in  donor‐ acceptor molecules, exemplified by para‐nitroaniline.   We will  review the applications of ATDA to prediction of the nonlinear  optical properties in the next Subsection.    2.4. Application  of ATDA  to predict  two photon absorption  profiles for fluorene­based chromophores  2PA  is  an  electronic  excitation  process  involving  simultaneous  absorption  of  two  photons.  There  are  a  wide  range  of  2PA  applications,  such  as  three  dimensional  data  storage, photonic devices, lithographic micro‐fabrication [114],  quantum  information  technology  [115],  optical  limiting,  two‐ photon pumped  lasing  in organic chromophores and quantum  dots  [115,116],  in‐vivo  bioimaging,  and  cell‐selective  photo‐ dynamic  therapy  [116].  Most  applications  require  chromophores with large 2PA cross‐sections to minimize laser  intensity  requirements  and  prevent  overheating  of  targets  [114].  To  design  more  efficient  2PA  chromophores,  it  is  important  to understand  their  structure/activity  relationships  (SARs).  Computer  modeling  of  2PA  spectra  facilitates  understanding  of  these  relationships  and  is  becoming  an  important part rational approach that may accelerate progress  in  the  chromophore design  [90].   Accurate predictions of 2PA  spectral  profiles  would  greatly  assist  in  the  design  of  more  effective 2PA chromophores while eliminating poor candidates  from  the  synthetic  pipeline.  The  goal  of  this  study  is  improvement  in  quantitative  predictions  of  2PA,  as  well  as  development of qualitative  tool  to  assist  in understanding  the  relations between the electronic structure of the chromophores  and 2PA profiles in terms of the transition dipole moments and  state energies.  In  the  past  decades  several  research  groups  had  made  a  strong effort aimed at the development of new compounds with  large  2PA  cross‐sections.  The  main  guiding  principle  used  in  those studies  involved the electron transfer between electron‐ donor  (D)  and  electron‐acceptor  (A)  moieties  attached  symmetrically  or  asymmetrically  to  the  π‐conjugated  bridge.  The  fluorene  fragment  in  particular  was  found  to  be  a  good  example  of  the  π‐conjugated  bridge  due  to  the  highly  delocalized  π‐system  delocalized  over  the  two  benzene  rings  held together at a nearly coplanar orientation by the methylene  bridge [117]. D‐π‐A, D‐π‐D, or A‐π‐A molecular structures have  been  proposed  and  studied  both  theoretically  and  experimentally. In studies by Belfield et al. fluorene derivatives  have been extended to D‐π‐π‐A and A‐π‐π‐A types with the aim  of  increasing  2PA  absorption  cross‐sections  [118‐123].  However, the choice of functional groups and linkages, the most  appropriate  for  developing  chromophores  with  large  2PA  characteristics, is still under active investigation [124‐135].  In  order  to  accelerate  experimental  efforts  based  on  the  traditional  trial  and  error  approach,  a  quantitative  understanding  of  the  trends  in  dependence  of  the  2PA  cross‐ section  on  a  molecular  structure  would  be  clearly  beneficial.  Two major approaches have been used to accomplish that goal.  The  first  one  is  based  on  essential  state  models  (three‐state,  four‐state, etc.). Parameters of these models (such as excitation  energies  and  transition  dipoles)  are  adjusted  to  fit  experimental data. These parameters are then correlated with  details of the molecular structure (π‐conjugated chain lengths,  donor/acceptor  strengths,  etc.).  Another  approach  consists  of  quantitative prediction of 2PA cross‐sections at a chosen  level  of theory, followed by analysis of the physical principles of the  major  contributions  into  this  property.  The  levels  of  theory,  used for 2PA predictions cover a wide range.  In  recent  quantum  chemical  calculations  performed  on  conjugated  chromophores  have  shown  that  a  substantial  symmetric  charge  redistribution  that  occurs  upon  excitation  may account  for a heightened sensitivity  to 2PA events  [118].  In  their  work  Bredas  et  al.  established  good  agreement  between the peak values of 2PA cross‐sections measured with  femto‐second  pulses  and  those  calculated  with  the  semi‐ empirical  intermediate  neglect  of  differential  overlap  Hamiltonian  with  multi‐reference  double‐configuration  interaction (INDO‐MRD‐CI) scheme based methods. Aside from  the donor‐acceptor configuration of 2PA active chromophores,  it  was  also  established  by  Bredas  et  al.  that  increasing  the  length  and  charge  transfer  of  the  molecules  results  in  an  increase  in  2PA  cross‐sections  and  may  also  result  in  a  significant shift of 2PA to longer wavelengths [118]. In another   150  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161      Table 2. Excitation energies ∆E (eV), permanent dipole moments, and leading electronic configurations for the lowest excited states of para‐nitroaniline at the  EOM‐CCSD theory level and four formalisms at the TD‐DFT level. The % of the double excited character is reported for EOM‐CCSD states and the numbers of the  excited state in each symmetry species are reported for TD‐DFT calculations. The ground state dipole moment is calculated to be 2.66 and 2.86 a.u. at CCSD and  DFT levels of theory, respectively.  Symmetry  EOM­CCSD  ∆E, eV  Configurations  %D    μ, a.u.    A1 (1)  4.92  (1(3)‐1'(3))+(2(4)2(4)‐1'(3)1'(3))‐(1(3)1(3)‐1'(3)1'(3))‐(1(3)1(4)‐1'(4)1'(3)) 0.9  5.65 6.99  (1(3)‐2'(3))+(1(4)‐1'(4))‐(1(3)1(3)‐1'(3)1'(3)) 0.3  4.64 7.48  (1(4)‐1'(4))‐(1(3)1(3)‐1'(4)1'(4))+(3(3)1(3)‐1'(4)1'(4)) 1.0  3.22 8.08  (2(3)‐1'(3))+(1(3)‐2'(3))‐(1(3)1(3)‐1'(3)1'(3)) 2.3  4.60 8.42  (2(4)‐1'(4))+(1(3)2(4)‐1'(3)1'(4)) 1.5  ‐1.59 B2 (2)  4.74  (1(3)‐1'(4))+(1(4)‐1'(3))‐(1(3)1(3)‐3'(3)1'(4)) 0.3  3.25 6.31  (1(4)‐1'(3))‐(2(4)‐1'(3))‐(1(3)‐1'(4))+(1(3)2(4)‐1'(3)1'(3)) 0.3  4.03 6.81  (2(4)‐1'(3))+(1(4)‐1'(3))‐(1(3)2(4)‐1'(3)1'(3)) 1.6  2.49 8.33  (2(3)‐1'(4))+(1(4)‐2'(3))‐(1(3)1(4)‐1'(4)1'(4)) 0.3  5.04 8.83  (1(4)‐2'(3))‐(2(3)‐1'(4))+(1(3)1(4)‐1'(4)1'(4)) 1.4  2.88 B1 (3)  4.76  (1(1)‐1'(3))+(1(1)‐2'(3))‐(1(2)2(4)‐1'(3)1'(3))+(1(2)2(4)‐1'(3)2'(3)) 0.9  2.03 6.58  (1(3)‐1'(1))+(2(3)‐1'(1)) 0.0  ‐0.16 7.69  (1(3)‐2'(1))  0.0  3.02 8.04  (1(2)‐1'(4))+(1(2)1(3)‐1'(4)1'(3)) 1.6  ‐1.54 8.56  (2(2)‐1'(4))  0.0  2.86 A2 (4)  4.26  (1(2)‐1'(3))+(1(2)‐2'(3))‐(1(1)2(4)‐1'(3)1'(3)) 0.5  1.92 7.61  (1(3)‐1'(2))  0.0  0.63 8.18  (1(1)‐1'(4))+(1(1)1(3)‐1'(4)1'(3)) 0.6  ‐0.75 8.33  (2(2)‐1'(3))  0.0  3.52 8.64  (1(4)‐2'(1))  0.0  1.19   Symmetry  TD­DFT  Root  ∆E, eV  Configurations  μ, a.u.  Numeric  Unre­laxed  ATDA  QR­DFT  A1 (1)  1  4.09 (1(3)‐1'(3))  4.93 7.61  5.56  4.93 2  6.21 (1(3)‐2'(3))‐(1(4)‐1'(4)) 3.59 4.42  3.66  3.60 5  6.98 (1(4)‐1'(4))+(1(3)‐2'(3))+(2(4)‐1'(4))+(2(3)‐1'(3)) 3.05 3.60  3.33  3.05 4  6.62 (2(3)‐1'(3))  5.06 5.96  4.34  4.61 3  6.52 (2(4)‐1'(4))  ‐0.30 ‐2.42  0.24  0.14 B2 (2)  1  4.59 (1(3)‐1'(4))+(1(4)‐1'(3)) 3.72 4.68  3.82  3.72 2  5.12 (1(4)‐1'(3))‐(1(3)‐1'(4)) 4.74 6.41  4.69  4.74 3  6.27 (2(4)‐1'(3))+(2(4)‐2'(3)) 1.03 0.97  0.79  1.03 4  7.23 (1(4)‐2'(3))‐(2(3)‐1'(4)) 3.69 4.58  3.74  3.69 6  8.00 (2(3)‐1'(4))+(2(4)‐2'(3))+(1(4)‐2'(3)) 2.12 1.68  2.18  2.39 B1 (3)  1  4.44 (1(1)‐1'(3))  2.11 1.50  2.16  2.11 2  5.96 (1(3)‐1'(1))  ‐0.40 ‐1.15  0.83  ‐0.40 4  6.92 (1(3)‐2'(1))  2.77 3.26  3.06  2.77 3  5.99 (1(2)‐1'(4))  ‐0.15 ‐3.27  ‐0.21  ‐0.15 6  7.91 (2(2)‐1'(4))  2.83 2.86  2.86  2.83 A2 (4)  1  3.96 (1(2)‐1'(3))  1.99 1.27  2.04  1.99 3  6.89 (1(3)‐1'(2))  0.55 0.19  1.52  0.53 2  6.35 (1(1)‐1'(4))  ‐0.01 ‐2.91  ‐0.02  ‐0.01 5  7.03 (2(2)‐1'(3))  4.43 5.83  4.35  4.45 7  8.09 (1(4)‐2'(1))  1.87 1.61  ‐0.63  1.88 rmsd    from QR‐DFT  0.60 2.08  0.70  ‐   from EOM‐CCSD  0.98 1.58  0.96  0.90     study  by  Agren  et  al.  [117]  four  lowest  excited  states    of  π  conjugated systems were theoretically investigated using using  the ab  initio response  theory. The materials were  synthesized  and characterized by Kim et al. [123] and Ventelon et al. [122]  Agren   et  al.   showed   that   their   theoretical   findings   were  consistent  with  the  correlation  between  large  2PA  cross‐ sections  and  a  π  center.  Although  1PA  spectra  were  strongly  correlated to the molecular length, this was not always the case  for  2PA  in  the  visible  domain  [117].  Yet  another  study  by  Fabian  et al.  established  that  spectral  absorption  features  are  reasonably well reproduced by the approximate semi‐empirical  MO‐CI methods based on the NDO (ZINDO/S), however the TD‐ DFT  response  theory  proved  to  be  superior  over  semi‐ empirical methods  [136].  Since  then  Hales  et  al.  showed  that  2PA spectra for symmetric and asymmetric fluorene derivative  compounds  exhibit  an  intermediate  resonant  enhancement  of  nonlinearities,  with  an  order  of  magnitude  enhancement  for  asymmetric  cases,  when  compared  to  degenerate  2PA.  Semi‐ empirical methods  that we  based  on  on  INDO‐MRD‐CI  theory  level and implemented a simplified three level model, were also  shown  to  provide  additional  insight  into  the  mechanisms  governing  2PA  events  [137].  Several  groups  published works  investigating  the  structure‐activity  relationships  (SARs)  responsible for these 2PA characteristics.  In  this Subsection we review  the  recent prediction of 2PA  cross‐sections  published  recently  [138].  The  conjugated  chromophores selected  in  that study to  illustrate  the accuracy  of  ATDA  method,  are  presented  in  Scheme  3.  Theoretical  models  of  these  were  derived  by  truncation  of  the  aliphatic  chains and replacing them with to methyl groups in the original  experimental  structures.  The  abbreviations  of  the  model  molecules  and  the  systematic  names  of  the  corresponding  experimentally studied ones are: 1a: 2,7‐Bisbenzothiazolyl‐9,9‐ didecylfluorene;  1b:  2,7‐Bis[4‐(9,9‐didecylfluoren‐2‐yl)vinyl]‐ phenylbenzothiazole;  1c:  9,9‐didecyl‐2,7‐bis(N,N‐phenyl  amino)‐fluorene;  1d:  (7‐benzothiazol‐2‐yl‐9,9‐didecylfluoren‐ 2‐yl)diphenylamine;  1e:  (7‐[2‐(4‐Benzothiazol‐2ylphenyl)  vinyl]‐9,9  didecylfluorene‐2yl)diphenylamine.  These  comp‐  ounds  were  experimentally  synthesized  and  characterized  by  Belfield et al. as a model compounds for possible applications in  two  photon  microfabrication,  two  photon  photochemical  transformations, non‐destructive 3‐D multiphoton fluorescence  imaging,  and  photodynamic  therapy  [119,137,139].  These  compounds exhibit  large  (600GM)  cross‐sections  for 1a while  studying  the  design  of  rigid‐rod  polymers  while  2PA  cross‐ Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  151  sections  of  the  polymers  were  reduced  by  aggregation  [140].  Compound  1a  additionally  exhibited  a  large  fluorescence  quantum yield. The good chemical, thermal, and photochemical  stability,  combined  with  desirable  one‐  and  two‐photon  absorption  and  luminescence  properties,  stand  out  as  characteristics of this chromophore as a promising material for  two‐photon based technologies [141]. Compound 1d has been  previously  implemented  for  in  vivo  2PA  biomedical  imaging  applications,  as  a  fluorophore  dye  used  for  staining  rat  cardiomyoblast  cells  (H9c2),  due  to  its  high  photostability,  fluorescence quantum yield, and two‐photon absorption cross‐ section  over  the  tunable  range  of  commercially  available  Ti:sapphire lasers [142]. Additionally, 1d has been investigated  as  a  potential  2PA  free‐radical  photo‐initiator  for  three‐ dimensionally  resolved  polymerization,  resulting  in  intricate  micro‐fabrication and imaging [143].  Of  the  compounds  in  Scheme  3,  1c  and  1d  have  been  previously  studied  theoretically  by  Hales  et  al  using  INDO‐ MRD‐CI  with  geometries  optimized  with  the  AM1  semi‐ empirical  Hamiltonian.  All  their  calculations were  carried  out  on  isolated  molecules  and  showed  strong  qualitative  and  quantitative agreement with experimentally generated spectra  [137]. The three‐level model developed in their study provided  further  insight  into  the  mechanisms  governing  the  nonlinear  activity  of  2PA  active  chromophores  in  relation  to  the  description of the molecular states. A second study carried out  by  Day  et  al.  implemented  linear  and  quadratic  response  density  functional  theories  to  calculate  the  photo‐physical  properties of D‐π‐A molecules  including 1d. Their comparison  of the two‐state approximation and with calculation of 2PA via  the  SOS  method  including  higher  energy  states  drew  a  conclusion  that  the  inclusion  of  higher  energy  states  was  necessary in description of 2PA [144,145].  In  this  Subsection  we  also  use  SOS  approach,  calculate  state‐to‐state  transition  dipole  moments i nm using  ATDA  as  detailed  in  Subsection  2.1  and  employ  them  to  identify  the  essential  states  governing  the  2PA  process.  Our  aim  is  to  validate ATDA by using  these approximate i nm   to predict  the  resonant  2PA  cross‐sections  with  SOS  models  and  compare  them to CEO results as well as experimental values.  Linear  response  TD‐DFT  allows  one  to  obtain  transition  dipole moments and transition energies from the ground state  to  a manifold  of  excited  states.  These  values  are  sufficient  to  calculate  a  1PA  spectrum.  Predictions  of  2PA  are  more  involved.  For  a  single  linear  polarized  laser  beam,  the  2PA  cross‐section,  averaged  over  isotropic  molecular  orientation  can be expressed in cgs units as follows [146]:         3 , ** 2 23 2 22 15 8 ba ababbbaaPA MMMMg c              (16)    Here a,b = x,y,z, and M  is so‐called 2PA transition moment,  expressed as the SOS via transition dipoles between the excited  states:                N k kg i kg j fk j kg i fk ij i M    2 1                                            (17)    Here g(2ω) is the Lorentzian lineshape. In this expression,  the  summation  runs  over  all  the  excited  states;  and  the  transition  dipole  moments  between  the  f‐th  and  k‐th  excited  states (permanent dipoles, when f  = k) are    fkkfkf   00                                      (18)  To  predict  a  2PA  spectrum  one  needs  transition  dipole  moments  between  excited  states  and  permanent  dipole  moments for infinite number of excited states. These values can  be obtained only if one retains the terms quadratic in external  field  while  solving  the  time‐dependent  equations.  This  quadratic approximation uses the states obtained in the linear  response approximation as  the basis, and  involves summation  over  infinite number of  these  states. This  summation can also  be  replaced  by  linear  equations,  as  it  is  done  in  QR‐DFT.  Alternatively,  one  can  carry  out  explicit  summation  over  a  sufficient number of the lowest excited states. The state to state  transition dipoles and permanent excited state dipoles can also  be  approximated  using  ATDA,  and  substituted  in  eq.  (16)  to  predict  2PA.  Transition  energies  and  ground  to  excited  transition  dipoles  remain  in  ATDA  identical  to  the  ones  obtained  in  full  linear  response  TD‐DFT  calculations.  At  the  same  time,  ATDA  allows  calculations  of  the  state  to  state  transition dipoles without solving the equations of full QR‐DFT.  Thus,  it  can  be  used  for  a  computationally  inexpensive  prediction  of  linear  and  nonlinear  properties  of  organic  molecules.  This  ATDA/SOS  approach  opens  new  venues  for  interpretation  of  2PA  properties  in  terms  of  molecular  electronic structure and can be used for rational design of 2PA  chromophores.    1a       1b   1c      1d      1e    Scheme 3. Structural formulas of the molecules studied.    All molecular structures were optimized at the HF/STO‐3G  theory  level,  which  favored  planar  geometry  of  conjugated  molecules and was shown [147] to give the best agreement for  the bond lengths as compared to the results of X‐Ray diffraction  experiments  for  stylbene  and  its  three  derivatives.  The  optimized  geometries  were  confirmed  by  the  absence  of  imaginary  frequencies  in  the  following  normal  mode  calculations.  The  single  point  energy  and  transition  dipole  calculations were performed at the TD‐B3LYP/MIDI! level.  Transition density matrices for the 24 lowest excited states,  as well  as  Kohn‐Sham operators  on  these  transition  densities  were  printed  out.  Contributions  of  the  second  and  third  derivatives  of  the  exchange‐correlation  potentials  into  Kohn‐ Sham  operators,  and  operators  on  the  pair  combinations  of  transition densities were neglected. The commercially available  computational  program  Gaussian98  [113]  was  modified  as  described  in  previous  studies  [148]  in  order  to  enable  this  printout.  The  frequency‐dependent  orientationally  averaged  first‐  and  third‐order  polarizability  tensors  were  generated  from the generated matrices using eq. (16), as implemented in  CEO program [74]. Habitual empirical  linewidth of 0.1 eV was  152  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161    used for both 1PA and 2PA. To analyze the electronic structure  of the excited states we used natural transition orbitals (NTO),  which  diagonalize  the  transition  density  matrix,  and  give  the  best  representation  of  the  electronic  excitation  in  single‐ particle  terms  [149].  The  graphical  software  XCrysDen  [150]  was used to plot NTOs.  We  present  2PA  resonant  energies  and  cross‐sections  in  Table 3 [138]. Figure 3 represents canonical orbitals, transition  dipoles and excitation energies for compound 1a. For two of the  molecules,  the  1PA  and  2PA  profiles  obtained  with  both  SOS  and  CEO  formalisms  are  presented  in  Figure  4,  along  with  linear spectra, and the natural transition orbitals. SOS and CEO  results (marked by solid lines on 2PA spectra in Figure 4) are in  excellent  quantitative  agreement  with  each  other,  which  provides a validation for the ATDA/SOS method. Predicted 2PA  profiles  also  agree well with  experimental  ones. Experimental  measurements  of  their  spectral  properties  were  reported  in  [140‐142,151]. While  for most molecules  agreement  between  resonant  maxima  is  better  than    0.1  eV,  in  the  case  of  1e  theoretical bands are red‐shifted relative to experimental ones  by approximately 0.44 eV (55 nm).  This can be explained by a  greater  conformational  flexibility  of  the  longer  conjugated  chain,  and  by  a  blue  shift  of  the  absorption  spectra  for  non‐ planar  conformations.  Overall,  agreement  with  experiment  validates  the  use  of  TD‐DFT  as  a  part  of  rational  design  strategies  directed  toward  new  and  improved  two‐photon  absorbing materials.        Figure 3. Canonical  orbitals  and  orbital  configurations  of  2PA  important  excited states, their excitation wavelength (), and transition dipoles ().    We  reviewed  theoretical  study  [138]  of  five  conjugated  chromophores experimentally shown to have large two photon  absorption cross‐sections. The use of the third order response  formalism within TD‐DFT to calculate the frequency‐dependent  second order hyperpolarizability in both SOS and CEO methods  to  describe  2PA  cross‐sections  resulted  in  very  similar  numerical  values.  The  interpretive  capabilities  of  these  methods  were  different,  however.  While  CEO  expressions  do  not lend themselves easily to a quantitative analysis, SOS ones  can  be  simplified  to  essential  state  models  and  employed  to  identify  2PA  resonant  states  and  interpret  the  relationships  between  the  electronic  structure  and  2PA  profiles.  Natural  Transition orbitals were also helpful to compare the electronic  structure of  the  linear and  two‐photon absorbing  states.  State  to  state  transition  dipole  moments,  necessary  for  SOS  expressions  are  calculated  with  ATDA  and  used  to  describe  two‐photon  processes.  Numerical  values  of  the  cross‐sections  obtained in SOS and CEO were found to be in good quantitative  agreement  with  each  other.  Both  SOS  and  CEO  results  are  in  good agreement with experiment. This validates the use of TD‐ DFT as a part of rational design strategies directed toward new  and  improved  two‐photon absorbing materials  for bioimaging  and optical data storage.    3. Density functional theory approach to photochemistry    3.1. Theory    In  this  Section we  review DFT‐based methods  that  can be  used  to  simulate  photochemical  processes,  such  as  light  initiated  bond  breaking  and  formation.  This  task  is  far  from  trivial,  as  KS  formalism  of  DFT  was  developed  for  non‐ degenerate  cases;  it  breaks  down  for  systems  with  strong  diradical  character  and  degeneracy  of  the  electronic  levels,  such as CIX geometries. However, the static (also known as left‐ right)  electron correlation  can be  taken  into  account by using  different  orbitals  for  different  spin.  This  approach,  known  as  unrestricted  KS  formalism  is  known  to  yield  a  qualitatively  correct description of the bond breaking [153]. Alternatives to  the  broken  spin‐symmetry  KS  solution  exist,  including  the  multireference DFT  [154],  the  local HF method  [155],  and  the  ensemble KS approach [156]. These new approaches, however,  require  development  of  new  exchange‐correlation  functionals  to obtain accurate numerical results.  Excited  states,  on  the  other  hand,  require  the  restricted  formalism to avoid a heavy mixing of the higher spin states into  description  of  the  excited  singlet.  Although  TD‐DFT  was  suggested  on  the  unrestricted  KS  reference  [157],  this  is  considered  to  be  incorrect  in  a  rigorous  theory  [158].  One  possible  approach  to  PES  of  the  excited  states  can  be  formulated  by  adding  the  excitation  energies  obtained  in  the  RTD‐DFT  formalism  to  the  ground  state  energies  calculated  with  the  unrestricted  KS method.  Thereafter  we  will  refer  to  this  approach  as  to  RTD‐UDFT.  Although  for  the  photoswitching  systems  considered  in  this  review  the  difference in the ground state energies obtained with restricted  and  unrestricted  KS  formalisms  are  close  to  within  20  kcal/mole or less, we will numerically show that this difference  is  sufficient  to bring excited states PES  to agreement with  the  results obtained at the higher theory level, when available.  At the first glance, the combination RTD‐UDFT may appear  contradictory.  Let  us  remind  the  reader  that  the  KS  approximation  to  the  DFT  is  semiempirical  by  virtue  of  the  exchange‐correlation  potential  used,  including  frequency  (in)  dependence  of  this  functional  to  describe  the  excited  states.  The  broken  symmetry  KS  formalism  compensates  the  deficiency of electron correlation description at the expense of  spin contamination. The symmetry dilemma (broken symmetry  or  restricted)  presents  the  choice  of  a  lesser  evil,  and  this  choice need not be the same for the ground and excited states.  While  spin  contamination  is  moderate  and  tolerable  for  the  ground state, we found it to be prohibitive for the excited state  perhaps  due  to  the  close  proximity  of  the  triplets  and  uncontrolled spin contamination.  Another theoretical development, necessary to describe the  region of CIX is related to the double excited states, missing in  the  adiabatic  linear  response  TD‐DFT  approximation  [47].  Mixing of  the double excited states  to  the  linear response TD‐ DFT states through the so‐called “dressed TD” was introduced  recently  by  Burke,  Maitra  et  al.  [159]  as  a  result  of  the  particular  frequency dependence  for  the  exchange‐correlation  potential. However, a much more conceptually simple solution  exists  in  the  CEO  formalism  [73‐75],  where  doubly  excited  states  appear  in  the  second  order  as  simple  products  of  the  excitations obtained at the linear response level. This fact was  used to propose ATDA, which we discussed in details in Section  2.  State‐to‐state  transition  dipoles  for  butadiene  and  permanent  dipoles  for  para‐nitroaniline  obtained  in  ATDA  were  compared  with  their  values  calculated  with  more  sophisticated  multireference  methods.  We  also  demonstrated  good  accuracy  of  the  ATDA,  comparing  one‐  and  two‐photon  absorption  spectra  for  several  chromophores  with  experimental data. Below we will review the studies  [30] and     Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  153        1a    1d        S1 particle    S1 hole    S4 particle    S4 hole    S1 particle    S1 hole    S3 particle    S3 hole  Figure 4.  Structural  formulas  (top  row),  1PA profiles  (row 2),  2PA profiles  (row 3),  and  isosurfaces  for natural  transition orbitals  (bottom  row)  for  studied  conjugated  chromophores.  Diamonds mark  the  experimental  profiles  taken  from  [140]  for 1a,  and  from  [137]  for 1d,  solid  lines  correspond  to  theoretical  predictions with SOS and CEO formalisms.    Table 3. Excitation energies (∆E, eV), wavelength (, nm), and 2PA cross‐sections  (, GM) for the electronically excited states in the molecules studied.  Experimental data are taken from (a) [140], (b) [137], (c) [152].   Compound  State  ∆Ecalc  ∆Eexp  λexp    Compound  State  ∆Ecalc  ∆Eexp  λexp  calc  exp    calc  exp   1a(a)  S1  3.59  3.41  364  1b(b)  S1  3.06  3.08  403 ‐  ‐  ‐  ‐ S2,4  4.28  4.27  290 S2  3.46  ‐  ‐ 324  437  711  ‐  1d(c)  S1  3.08  3.21  387 S4  3.99  4  310 65  73  ‐  ‐ S3  4.10  ‐  ‐ S11  4.45  ‐  ‐ 151  ‐  486  ‐ S4  4.12  ‐  ‐  1c(b)  S1  3.38  3.1  400 151  ‐  ‐  ‐  1e(c)  S1  2.83  3.24  383 S3  3.88  4  310 154  162  126  89 S2  3.56  4.03  306 S5,6  4.01  ‐  ‐ 133  ‐  ‐  ‐ S6,7,8  4.25  ‐  ‐ S15  4.6 ‐  ‐ 445  ‐  162  ‐       154  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161    [160], where  ATDA‐UDFT  produced  the  accurate  energies  for  the HOMO‐LUMO double excited state in the entire range of the  bond breaking reaction coordinate C1…C6 for some benchmark  molecular  photoswitches,  provided  that  the  molecular  geometry corresponds to that state.  Since analytical gradients in the ATDA‐UDFT approach are  not  yet  implemented  in  the  computer  codes,  the  studies  reviewd  below  use  an  approximation  called  Slater  transition  state (STS) method to optimize geometry of the excited states.  In  this  method,  half  an  electron  is  promoted  from  HOMO  to  LUMO  and  self‐consistency  is  achieved  with  these  fractional  orbital  occupations  [161].  STS  is  known  to  be  a  good  approximation  to  the  corresponding  ΔSCF  excitation  energy  [162,163].  It  can  be  further  extended  to  the  modified  linear  response  DFT  method  [164],  which  yields  a  considerable  improvement in description of the charge‐transfer and Rydberg  states,  compared  to  the  linear  response  TD‐DFT  approach.  A  practical  advantage  of  STS  is  an  easier  SCF  convergence,  compared  to  the  excited‐state  SCF  convergence,  which  often  presents a major problem [161] [23].  Studies  [30]  and  [160]  used  the  Gaussian  2003  Rev.  E1  suite  of  programs  [112],  and  hybrid  meta‐GGA  exchange‐ correlation functional M05‐2X from Truhlar’s group [165]. This  functional  was  designed  for  accurate  description  of  both  equilibrium  geometries  and  transition  states.  Relaxed  energy  scans  (RES)  were  built  by  means  of  the  constrained  optimization technique and used as approximations to the true  minimum energy pathways  [166].  The pericyclic  C1…C6 bond  (reaction  coordinate) was  fixed  in  the  range  from 1.4  to  3.6Å  with 0.1Å step size. The ground singlet state 1A, was optimized  at  the  UM05‐2X/6‐31G  level  of  theory,  while  single  1B  and  double  2A  excited  singlet  state  geometries were  optimized  in  the STS method. The STS method was implemented using equal  fractional  occupation  numbers  for  HOMO  and  LUMO  in  the  alpha‐set  only  [to  approximate  geometry  of  the  single  excited  state 1B,  IOp(5/75=1, 76=2)], and both alpha and beta sets  to  approximate  geometry  of  the  double  excited  state  2A  [IOp(5/75=1,  76=2,  77=1,  78=2)]. While  the  STS‐UM05‐2X/6‐ 31G  theory  level  gave  reasonable  estimates  for  molecular  geometry, the SCF energies obtained at this level were found to  have  an  insufficient  accuracy.  Instead,  the  excitation  energies  were taken from single‐point calculations at the corresponding  optimized geometry in ATDA for the lowest single excited state  1B  and  the  lowest  double  excited  state  2A.  The  excitation  energies thus obtained at the ATDA‐M05‐2X/6‐31G level, were  added  to  the  ground  state  energies  obtained  at  UM05‐2X/6‐ 31G  level  of  theory.  Resulting  ATDA‐UM05‐2X/6‐31G//STS‐ UM05‐2X/6‐31G energies were plotted  in  Figure 6–Figure 10.  The KS orbitals in Figure 11 were generated using VRML option  in the Molden graphical visualization package [167].   3.2. Benchmark study of hexatriene    Cyclohexadiene/hexatriene  conversion  is  the  simplest  example  of  an  electrocyclic  reaction.  The  dynamics  of  the  cycloreversion  in  the CHD/HT  system was  repeatedly  studied  with  the  time‐resolved  ultrafast  spectroscopy  techniques  [19,168].  The  results  are  summarized  in  Ref.  [17].  While  a  solvent  may  stabilize  the  ionic  1B  state  compared  to  the  covalent  2A  and  1A  states,  this  does  not  have  an  appreciable  effect  on  the  CHD/HT  conversion  rate,  while  it  slows  down  conformational  isomerization  of  HT  considerably  [169].  To  describe  the  detailed  mechanism  of  this  photoprocess,  it  is  convenient  to  consider  state  correlation  diagrams,  briefly  introduced below.  The  CHD/HT  system  presents  a  textbook  case  of  an  electrocyclic reaction, considered in the very first publication of  orbital  symmetry  conservation  rules  by  Woodward  and  Hoffman [170]. Their pioneering work pointed out that in order  for the HT→CHD transition to have a low activation barrier in  the  ground  state  (be  thermally  allowed),  the  HOMO  should  have bonding character with respect to the C1–C6 interaction,  closing  the  ring.  This  condition  is  satisfied when HOMO  lobes  localized  on  C1  and  C6  atoms  overlap  in  phase  during  disrotatory motion of the terminal CH2 groups. Thus, the HOMO  symmetry  controls  stereochemistry  of  the  thermal  reaction  product when substituents are present. Promotion of electron  to  LUMO  in  the  photoactivated  process  inverses  symmetry  of  the wave  function.  In  this  case  conrotatory motion  of  HT  (or  CHD  in  cycloreversion  reaction)  is  necessary  for  lobes  on  C1  and  C6  atoms  to  overlap  in‐phase  (Figure  5)  and  create  no  potential  barrier  along  the  reaction  pathway,  making  the  reaction photochemically allowed.  Longuet‐Higgins and Abrahamson reformulated the rule in  terms  of  state  correlation  diagrams  [171]  (Figure  5,  center).  These diagrams consider three electronic states of the HT/CHD  system, classified as symmetric (1A, 2A) or antisymmetric (1B)  with  respect  to  the  axis  C2,  conserved  along  the  conrotatory  cyclization pathway. State 2A is the doubly excited state of the  reactant  (2  electrons  on  LUMO,  0  electrons  on  HOMO).  It  correlates with the ground state of the product, which has the  same  symmetry  (1A),  while  the  singly  excited  state  of  the  reactant  (1B)  correlates  with  the  singly  excited  state  of  the  product (HOMO and LUMO occupied with 1 electron each). It is  important to note that selection rules allow 1A→1B transitions  in  the  1PA  regime,  and  the  1A→2A  transition  in  the  2PA  regime. Along the reaction pathway 1A raises, and 2A drops in  energy. When approaching degeneracy, 1A and 2A mix to form  an  avoided  crossing.  This  region  forms  a  potential  barrier  on  the ground state surface, and the so‐called pericyclic minimum  on  the  doubly  excited  state  (2A)  surface.  Even  though  energy  diagrams  are  drawn  along  the  C2‐symmetric  coordinate  and  clearly  present  an  oversimplification,  qualitative  diagrams  based  on  semiempirical Molecular  Orbital  calculations  served  as  a  useful  guidance  for  the  rational  design  of  thermally  irreversible photochromic molecules [5].  It is worth noting that σ and σ* orbitals of CHD, involved in  the bond breaking process (Figure 5), are not frontier orbitals  (HOMO  and  LUMO),  unless  the  bond  is  elongated  enough.  As  was  stressed  by  Michl  [172],  the  photochemical  transformations  are  controlled  by  essential  orbitals,  which  often do not coincide with frontier orbitals. This means that the  photochemically excited  state  (usually  the  lowest  singlet  state  of essentially HOMO→LUMO character) may increase in energy  along the reaction pathway before  it can proceed down to the  pericyclic  minimum.  This  may  lead  to  an  activation  barrier  which will  slow  down  or  prohibit  the  photochemical  process.  Thus,  even  simple  orbital  considerations  can  be  used  for  semiquantitative prediction of the photochemical kinetics.  More accurate quantitative calculations, on the other hand,  involve plotting realistic potential energy surfaces (PES), which  until recently required the use of ab initio multireference‐based  quantum chemistry methods. Pioneering CAS study of HT/CHD  system was published by Robb, Olivucci et al.  [173‐175]. They  found  that  2A  and  1A  surfaces  touch  at  a  molecular  conformation of  tetraradical  character,  located  away  from  the  C2‐symmetric  coordinate,  and  including  partial  bonds  C1…C6  and C2…C6. They also found that account for dynamic electron  correlation is essential to correctly predict relative energies of  1B and 2A states.  Despite  the  fact  that  2A/1A  and  1B/2A  conical  intersections  complicate  the  energy  landscape  by  adding  an  extra  dimension,  qualitative  interpretation  given  by  the  state  correlation diagram along the symmetric coordinate still holds.  It  was  further  confirmed  at  the  high  level  (CAS  second  order  perturbation  theory  and  multireference  configuration  interaction)  by  building  ab  initio  PES  and  performing  two‐ dimensional  quantum  dynamics  [176,177]  on  this  PES.  A  barrierless  descent  on  the  excited  state  PES  was  found  to  determine ultrafast photoconversion between CHD and HT, and  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  155  the quantum yield of this process was primarily determined by  the location of 2A/1A CIX.      Figure  5.  Orbital  (top)  and  state  (center)  correlation  diagrams  for  photointiated  electrocyclic  conversion  of  CHD  to  HT  (bottom).  Lines  connect the orbitals as they continuously transform along the C2‐symmetric  reaction path. Different orbital occupations give rise to different electronic  states. The doubly excited state 2A of CHD correlates with the ground state  of HT and vice versa; when their energies approach each pother, the states  strongly mix  and  avoid  crossing;  a  true CIX occurs  away  from  symmetric  coordinate.    Relaxed  potential  energy  scans  obtained  in  this  study  are  plotted in Figure 6. They are in good agreement with high level  ab  initio  results.  As  the  C1…C6  reaction  coordinate  contracts  from  non‐bonding  distance  to  the  normal  covalent  bond,  the  bright 1B state (characterized by a large transition dipole from  the ground state),  is being monotonically stabilized  in energy,  starting  descent  to  the  pericyclic minimum.  At  the  same  time  the dark 2A state (with a negligibly small transition dipole from  the ground state) is being stabilized even faster, crosses the 1B  state surface, and forms the bottom of the pericyclic minimum.  It  is  worth  noting  that  the  geometry  of  the  ground  state  does  not  approximate  the  excited  state  geometry  accurately  enough  to  produce  a  reasonable  RES.  The  vertical  excitation  energy  curve,  plotted  in  Figure  6  and  representing  energy  of  the  1B  excited  state  at  the  ground  state  geometry,  displays  a  maximum  instead  of  the  pericyclic minimum,  and  contradicts  both  high  level  ab  initio  data,  and  the more  accurate  relaxed  ATDA‐UDFT data.  Qualitative  comparison between our modified DFT  results  and  a  state‐of‐the‐art  wavefunction  method,  multi  reference  second order perturbation theory, calculations of  the CHD/HT  system  is  presented  in  Figure  7.  Four  important  points  were  considered:  equilibrium  geometries  for  closed  and  open  isomers  (CHD  and  HT),  corresponding  to  the  Franck‐Condon  geometries of the ground states; and excited states 1B and 2A,  optimized into the respective pericyclic minima. As one can see,  ATDA‐UDFT  //STS‐DFT  at  the  M05‐2X/6‐31G  theory  level,  adopted  in  this  work,  almost  uniformly  overstabilizes  both  excited states by 0.4–1.0 eV, but retains the correct order of the  excited states, as compared to the multireference perturbation  theory results. To the best of our knowledge,  the correct state  ordering  in  this  system was  not  predicted  by  any  other DFT‐ based approach.      Figure 6. Relaxed potential energy scans for the ground (S0, 1A – solid line),  single excited (S1, 1B – long dashes) and double excited states (S2, 2A – short  dashes)  along  the  reaction  pathway  of  the  ring‐closing  C1…C6  bond  in  the  CHD/HT  system,  predicted  at  the  ATDA‐UDFT/6‐31G//STS‐DFT/6‐31G  theory  level, using  the M05‐2X exchange‐correlation  functional. The absence  of an appreciable energy barrier on the pathway from the 1B state of CHD in  the  Franck‐Condon  geometry  (left)  to  the  minimum  on  the  2A  surface  is  consistent with an ultrafast rate of the photoinitiated cycloreversion reaction  CHD→HT.  The  upper  dash‐dotted  line  corresponds  to  the  vertical  single  excitation at the optimized ground state geometry (S1, 1A).      Figure 7. Relative state energies (in eV) for the ground and two lowest singlet  excited  states  in  the  CHD/HT  system,  obtained  at  the  ATDA‐UDFT/6‐ 31G//STS‐DFT/6‐31G  theory  level,  using  the  M05‐2X  exchange‐correlation  functional.  The  results  of  high  level  ab  initio  calculations  from  Ref.  [176]  (multi‐configurational quasi‐degenerated perturbation theory) are shown for  comparison in parentheses.    3.3 Benchmark study of dithienylethene    Dithienylethene  (Scheme  1  with  X=S)  is  the  simplest  homolog of diarilethens,  an  important  class of  compounds  for  photoswitching  applications.  In  the  closed  conformation  its  thiophene rings are  joined through the central bond to  form a  conjugated  π‐electron  system.  In  the  open  form  the  π‐ conjugation  between  the  rings  is  disrupted  by  the  large  156  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161    dihedral  twists.  The  photochemistry  of  dithienylethene  upon  one‐photon  irradiation  was  studied  in  Ref.  [9]  using  a  molecular  mechanics–valence  bond  computation.  Dynamics  was  simulated  on  a  hydrocarbon  model  system  having  the  same potential energy surface topology as dithienylethene. The  critical  points  were  obtained  at  the  CAS(10,10)  theory  level  with  basis  sets  6‐31G  and 6‐31G*. Dynamics was  found  to  be  essential  to  explain  the  observed  photochemical  behavior  qualitatively. The authors claim, that the reaction path alone in  insufficient for that purpose, as initial relaxation for this system  was  not  directed  toward  CIX.  However,  the  shape  of  the  first  excited state RES was found to be useful in interpretation of the  ultrafast  ring‐closure  and  a  relatively  low  quantum  yield  observed  for  the  ring‐opening  reaction.  Stepwise multiphoton  processes  in  dithienylethene  were  studied  in  Ref.  [178]  The  RES of  several  lowest singlet excited states were calculated at  CAS(10,10)/6‐31G  and  CAS(12,14)/6‐31G  theory  levels.  In  all  three of the excited states reported, the potential surfaces had a  barrier  to  cycloreversion,  located  between  the  closed  ring  conformation  and  the  pericyclic  minimum.  Although  the  crossover between  excited  states was not  found  in  that  study  (presumably,  due  to  neglect  of  dynamic  electron  correlation),  the height of the barrier was found to be much lower in higher  lying  states.  This  fact was  used  to  explain  the  experimentally  observed  sequential  two‐photon  photoinitiation  (via  excited  state absorption) in this system.  Our  approximations of  the minimum energy pathways  for  the  ground  and  lowest  single  and  double  excited  states  are  plotted in Figure 8. One can see that RES of the 1B state has a  minimum  in  the  closed  form,  and  monotonically  rises  to  the  Franck‐Condon region of the open form. The doubly excited 2A  state, on the other hand, forms pericyclic minimum and crosses  below  the  1B  state  in  the  vicinity  of  CIX.  Therefore,  the  excitation  of  the  open  form  is  followed  by  the  ultrafast  barrierless  relaxation  into  the  pericyclic  minimum  along  1B  and  then  2A RES, while  the  excitation  of  the  closed  form will  populate  the  potential minimum  on  the  excited  state  surface.  The  conversion  of  the  excited  closed  form  into  the  pericyclic  minimum  must  first  overcome  a  small  (c.a.  5  kcal/mol)  potential  energy  barrier,  which  leads  to  relatively  slow  cycloreversion. The excited state absorption will then bring the  system  from  1B  to  2A  state,  followed  by  the  barrierless  relaxation  toward  the  CIX.  Thus,  RES  explain  both  slow  cycloreversion  and  ultrafast  photoswitching  upon  the  sequential two‐photon absorption.      Figure 8. Relaxed energy scan for the ground (S0, 1A), single excited (S1, 1B)  and  double  excited  states  (S2,  2A)  along  the  reaction  pathway  of  the  ring‐ closing C1…C6 bond  in  the dithioarylethene  system, predicted  at  the ATDA‐ UDFT/6‐31G//STS‐DFT/6‐31G  theory  level,  using  the  M05‐2X  exchange‐ correlation functional. The legend is the same as in Figure 6.        3.4. Application to the prototype photoswitching molecule    Compound  2b  (Scheme  4)  was  recently  synthesized  by  Belfield  et  al.  [152].  They  covalently  attached  the  diphenylaminofluorene moiety  to  the  compound 2a.  RESs  for  2a are plotted in Figure 9, and indicate that the excited state of  the open form relaxes barrierlessly to the pericyclic minimum,  while  the  closed  form  in  its  lowest  excited  state  is  separated  from  the  pericyclic  minimum  by  a  sizable  potential  energy  barrier. Therefore, ATDA‐UDFT//STS‐DFT method predicts 2a  to  undergo  ultrafast  photochemical  cyclization,  while  the  cycloreversion  is  expected  to  be  slow,  in  agreement  with  experimental data [179]. In the experiment, compound 2a was  doped  into  the  polymer  matrices  (polyvinylcarbazol  and  polymethylmethacrylate). In either case the rate constant of the  coloring  process  (cyclization)  was  found  more  than  twice  higher than the rate of the bleaching process (cycloreversion).  The absorption spectrum of diarylethene 2a in the closed form  has  two  bands  (at  470–550  nm  and  360–390  nm).  TD‐DFT  calculations for gas phase yield the absorption maximum with a  considerable blue shift, when compared to the experiment (the  lowest vertical singlet excitations are at 423 and 324 nm for the  closed ring geometry). Upon irradiation with visible light in the  experiment,  the  first  band  disappears  and  the  second  one  is  considerably  reduced  in  intensity.  This  is  reproduced  by  TD‐ DFT  calculations.  The  first  singlet  excited  state  for  the  open  form  2a  is  at  386  nm,  but  the  first  excitation  with  sufficient  oscillator strength is at 216 nm.   The  diphenylaminofluorene  substituent  was  designed  similar  to  other  fluorene‐based  donor‐fluorene‐donor  and  donor‐fluorene‐acceptor  chromophores,  previously  shown  to  have a large 2PA cross‐section [180,181]. Its attachment to the  photoswitching  core 2a  (Scheme  4) was  expected  to  produce  an  enhanced  photochromic  compound  designed  to  work  as  Write  Once,  Read  Many  (WORM)  photoswitch  in  the  two‐ photon absorbing regime.    S S N O O H    (2a)  S S N O O N   (2b)    S S N O O N F F F F F F F F F F FF F F F F    (2c)   Scheme 4. Open forms of considered compounds.     Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  157    Figure 9. Relaxed energy scan for the ground state and the lowest single and  double excited states in compound 2a (Scheme 4). The legend is the same as in  Figure 6. There is no barrier along the single excited state from the open‐ring  Franck‐Condon  region  toward  the  pericyclic  minimum.  Experimentally,  this  compound is photoreactive in the both closed and open forms, though the rate  of  the  coloring  process  (cyclization)  is  higher  than  the  rate  of  bleaching  (cycloreversion).    While  the  compound  2a  was  experimentally  found  to  be  photoreactive,  the molecule 2b  did  not  display  photochromic  properties.  ATDA‐UDFT//STS‐DFT  approach,  verified  in  the  Subsections  3.2  and  3.3,  was  applied  to  the  molecule  2b  in  order to predict the RESs. The resulting curves are reported in  Figure  10b.  As  one  can  see,  RES  of  the  lowest  excited  state  forms minima  for  the  both  open  and  closed  forms,  separated  from the pericyclic minimum by potential energy barriers of 4.7  and 8.0 kcal/mole respectively. As it was shown in the previous  Subsections,  the  barrier  for  the  closed  form  is  typical  for  the  other  diarylethene  derivatives.  However,  the  barrier  for  the  open  form  is  unexpected  and  is  responsible  for  the  loss  of  photoswitching properties.   In order to determine the nature of this barrier, KS orbitals  in  these  molecules  were  analyzed.  This  analysis  found  that  while  the  essential  orbitals  (antibonding  and  bonding  with  respect  to  forming the new C–C bond) are respectively HOMO  and  LUMO  in  molecule  2a,  they  are  HOMO–1  and  LUMO  for  molecule  2b.  At  the  same  time  HOMO  in  molecule  2b  is  primarily the lone pair of the nitrogen atom and is localized on  the  diphenylaminofluorene  pendant  group  (Figure  11),  designed to increase the 2PA cross‐section. In molecule 2a the  HOMO to LUMO excitation reduces the anti‐bonding effect and  creates a bonding one. This  increases the bond order between  carbon  atoms  closing  the  ring  from  –1  to  0  and  allows  the  descent  toward  the  pericyclic  minimum  on  the  excited  state  PES.  On  the  other  hand,  in molecule 2b  the  change  in  HOMO  occupation has no pronounced  effect  on  the  ring  closure.  The  antibonding electrons are kept on HOMO–1 upon the excitation.  The placement of an electron into bonding LUMO increases the  bond  order  from  –1  to  –1/2  with  net  antibonding  effect  remaining, which is apparently insufficient for initiation of the  ring closure reaction. Therefore,  the  lowest excited state  is no  longer photoreactive. Not  surprisingly,  its  energy  raises  along  the  reaction  coordinate  connecting  the  Franck‐Condon  region  and the pericyclic minimum. This is exactly what was observed  on RES plotted in Figure 10b.  Apparently,  in  order  to  restore  the  photoreactivity  in  the  prototype molecule one needs  to stabilize  the  lone pair of  the  donor  substituent.  This  can  be  accomplished  by  chemical  modifications  of  the  pendant  group  with  an  electronegative  substituent.  If  this  stabilization  is  sufficient,  the  donor  orbital  will have  the  lower energy  than both of  the essential orbitals,  and  they  become  the  HOMO‐LUMO  pair  once  again.  The  perfluorination  of  the  pendant  group  should  stabilize  all  the  orbitals,  localized  on  this  molecular  fragment.  Indeed,  it  was  shown  to  result  in  greated  stabilization  of  the  donor  substituent  lone pair,  relative  to  the pair  of  essential  orbitals.  The modified prototype  (molecule 2c)  is  shown on Scheme 4.  Comparison of the energy levels in the molecules 2b and 2c  is  presented  in  Figure  11. While  perfluorination  of  the  pendant  group  lowers  energy  for  all  the  orbitals  considered,  the  stabilization  is  the  greatest  for  the  lone  pair  orbital,  localized  primarily on this fragment.        (b)      (c)    Figure 10. Relaxed  energy  scan  for  photocyclization  in molecules 2b (b,  top) and 2c (c, bottom), according to our results at the ATDA‐UDFT//STS‐ DFT, M05‐2X/6‐31G level of theory. The legend is the same as in Figure 6.  While for 2b the both excited states exhibit potential energy barriers from  the  open  isomer  to  the  pericyclic  minimum  (b,  top),  there  is  no  such  a  barrier  in  case  of  2c  (c,  bottom).  This  indicates  that  an  ultrafast  photoconversion is expected for 2c.    The  plot  of  RESs  for  the  excited  states  in  molecule  2c  (Figure 10c) confirms the essential orbital analysis. The energy  of the lowest excited state is monotonically decreasing from the  Franck‐Condon  region  of  the  open  form  to  the  pericyclic  minimum with no barrier. Thus, the electronic structure of the  reactant can be used to predict the nature of the lowest excited  state of  the molecule, potential energy surface  for  this excited  state,  and  the  photochemical  kinetics  of  electrocyclization.  Based on the electronic structure of the systems considered in  this Subsection, a general design principle for the photochromic  molecules  was  formulated  [30]:  substituents  should  not  introduce  any  electron  energy  levels within  the HOMO‐LUMO  gap between the essential orbitals.  158  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161          Figure 11. Kohn‐Sham energy levels (a.u.) and orbital plots for open forms of the compounds 2b (left) and 2c (right) at the M05‐2X/6‐31G level of theory. While  the orbitals essential for the ring closure are HOMO–1 and LUMO in 2b, they are HOMO and LUMO in 2c. This is the reason why the lowest excited state of 2c is  photoreactive.        Figure 12. 1PA (bottom row) and 2PA (top row) spectra for molecules 2b (right column) and 2c (left column) calculated at the ATDA‐B3LYP/6‐31G//AM1 level  of theory. Vertical lines correspond to calculated energies of excited states.    The  theoretically  predicted  1PA  and  2PA  spectra  for  the  open  form  of  molecules  2b  and  2c  are  plotted  in  Figure  12.  Both molecules demonstrate moderate 2PA cross‐sections (ca.  30  GM)  in  the  range  of  excitation  energy  300  –  600  nm.  To  increase  the  2PA  cross‐section  some  additional  chemical  modifications  are  necessary.  This  can  be  a  subject  of  future  studies.    4. Conclusions    We  reviewed  the  current  state  of  DFT‐based  methods  developed  for  prediction  of  2PA  cross‐sections  and  potential  energy  surfaces  of  the  excited  states.  A  recently  proposed  ATDA‐DFT approach,  based on  the  approximate  second order  response  within  density  functional  theory  was  shown  to  be  capable  to  accomplish both goals. Apart  from  the quantitative  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  159  accuracy, this approach allows for qualitative interpretation of  the electronic structure that leads to a better understanding of  the  photophysical  and  photochemical  processes  involved.  ATDA‐DFT  approach  (a  posteriori  Tamm‐Dancoff  approximation) accounts for both single and double excitations  appearing in the first and second order time‐dependent DFT in  the  Coupled  Electronic  Oscillator  formalism.  The  permanent  and  transition  dipole  moments,  calculated  at  this  approximation,  reproduce well high‐level ab  initio values,  and  ATDA predicted 2PA cross‐sections are  in agreement with  the  third‐order  TD‐DFT  results.  Smooth  potential  energy  surfaces  were obtained by replacing the ground state energy component  of  the  excited  state with  the  one  obtained  in  the  unrestricted  broken  symmetry Kohn‐Sham  formalism  (termed here ATDA‐ UDFT). The importance of excited state geometry optimization  (as  opposed  to  the  habitual  use  of  unrelaxed  ground  state  geometries)  in  accurate  prediction  of  these  potential  energy  surfaces was confirmed. For the lack of analytical derivatives at  the  ATDA‐UDFT  theory  level  the  lowest  single  and  double  excited  state  geometry  is  approximated  using  the  Slater  Transition  State  method  (STS‐DFT).  The  combined  ATDA‐ UDFT//STS‐DFT  approach  was  shown  to  somewhat  underestimate  energies  of  the  both  excited  states  but  closely  reproduce  the  state  ordering  and  the  energy  crossovers  predicted  at  high‐level  multireference  perturbation  theory  level  for  the  prototypical  hexatriene/cyclohexadiene  system.  The  approach  was  also  able  to  explain  the  experimentally  observed  slow  photochemical  cycloreversion  rates  in  model  dithioarylethenes.  The  ATDA‐UDFT//STS‐DFT  approach  was  also  used  to  explain  the  loss  of  photochromic  reactivity  in  the  prototype  two‐photon absorbing photoswitch. The nature of the loss was  traced  to  the  presence  of  a  barrier  on  the  potential  energy  surface  of  the  lowest  excited  state.  The  Kohn‐Sham  orbital  analysis  revealed  this barrier  to be  the  result  of  the  lone pair  orbital  of  the donor  chromophore  group  appearing  inside  the  HOMO‐LUMO  gap  of  the  essential  orbitals  in  the  photoswitching molecular core, so that the lowest excited state  changed  from  the  photoreactive  type  to  the  charge  transfer  type. In order to restore the photoreactivity the stabilization of  the lone pair orbital to the energy below the HOMO‐LUMO gap  was  suggested by means of perfluorination of  the  two‐photon  absorbing chromophore pendant group. This  substitution was  found  to  restore  the  order  of  the  essential  orbitals  and  eliminate  the  energy  barrier  on  the  excited  state  potential  surface.  Based  on  these  findings,  the  design  principle  for  the  new  photochromic  molecules  was  formulated:  substituents  should  not  introduce  any  electron  energy  levels  within  the  HOMO‐LUMO  gap  between  essential  orbitals  of  the  photoreactive  fragment.  This  principle  will  assist  in  future  development  of  the  materials  for  three‐dimensional  optical  data  storage  devices with  unprecedented  storage  capacity  for  advanced applications in Information Technology.  The ATDA‐DFT  approach was  also  found  to  predict  state‐ to‐state transition dipoles and permanent excited state dipoles  accurately, as compared to the correlated wavefunction results.  The  two‐photon  absorption  cross  sections  predictions  were  also  accurate,  as  compared  to  the  experimental  data  for  a  benchmark  set.  Thus,  ATDA‐DFT  provides  a  tool  that may  be  used  in  design  of  more  efficient  two‐photon  absorbing  chromophores.  The  methods  reviewed  here  are  capable  to  accurately  describe  both  nonlinear  optical  and  photochemical  properties of the material; they call for broader application and  further development.    Acknowledgements    This  work  is  supported  in  part  by  the  National  Science  Foundation  (CCF  0740344).  Research  was  performed  in  part  using  (1)  Stokes HPCC  facility  at  UCF  Institute  for  Simulation  and Training  (IST),  (2) Bethe SMP server  at UCF NanoScience  Technology  Center  (NSTC),  and  (3)  the  National  Energy  Research Scientific Computing Center (NERSC) a DOE Office of  Science user  facility at Lawrence Berkely National Laboratory.  AEM  acknowledges  ACS  COMP  Hewlett‐Packard  Outstanding  Junior Faculty Award presented for this work at the Fall 2008  Annual Meeting of the American Chemical Society.     References  [1]. Kaneko, M. MRS Bull. 2006, 31, 314‐317.  [2]. Wang,  J.  G.;  Sun,  C.  J.;  Hashimoto,  Y.;  Kono,  J.;  Khodaparast,  G.  A.;  Cywinski, L.; Sham, L.  J.; Sanders, G. D.; Stanton, C.  J.; Munekata, H.  J.  Phys.: Condens. Matter 2006, 18, R501‐R530.  [3]. Zhou, G. F. Mater. Sci. Eng., A 2001, 304, 73‐80.  [4]. Kawata, S.; Kawata, Y. Chem. Rev. 2000, 100, 1777‐1788.  [5]. Irie, M. Chem. Rev. 2000, 100, 1685‐1716.  [6]. Parthenopoulos, D. A.; Rentzepis, P. M. Science 1989, 245, 843‐845.  [7]. Saita, S.; Yamaguchi, T.; Kawai, T.; Irie, M. Chem. Phys. Chem. 2005, 6,  2300‐2306.  [8]. Shipway, A. N.; Greenwald, M.;  Jaber, N.; Litwak, A. M.; Reisman, B.  J.  Jpn. J. Appl. Phys., Part 1 2006, 45, 1229‐1234.  [9]. Boggio‐Pasqua, M.; Ravaglia, M.; Bearpark, M. J.; Garavelli, M.; Robb, M.  A. J. Phys. Chem. A 2003, 107, 11139‐11152.  [10]. Koppel, H.; Doscher, M.; Mahapatra, S. Int. J. Quantum Chem 2000, 80,  942‐949.  [11]. Kendrick, B. K. J. Chem. Phys. 2001, 114, 8796‐8819.  [12]. Hack, M. D.; Jasper, A. W.; Volobuev, Y. L.; Schwenke, D. W.; Truhlar, D.  G. J. Phys. Chem. A 2000, 104, 217‐232.  [13]. Jasper, A. W.; Truhlar, D. G. J. Chem. Phys. 2005, 122, 044101‐16.  [14]. Vallet,  V.;  Lan,  Z.  G.;  Mahapatra,  S.;  Sobolewski,  A.  L.;  Domcke,  W.  Faraday Discuss. 2004, 127, 283‐293.  [15]. Garavelli, M. Theor. Chem. Acc. 2006, 116, 87‐105.  [16]. Truhlar, D. G.; Gordon, M. S. Science 1990, 249, 491‐498.  [17]. Fuss,  W.;  Hering,  P.;  Kompa,  K.  L.;  Lochbrunner,  S.;  Schikarski,  T.;  Schmid, W. E.;  Trushin,  S. A. Ber. Bunsen­Ges. Phys. Chem 1997, 101,  500‐509.  [18]. Bernardi, F.; Olivucci, M.; Robb, M. A. Chem. Soc. Rev. 1996, 25, 321‐ 328.  [19]. Fuss, W.; Lochbrunner, S.; Muller, A. M.; Schikarski, T.; Schmid, W. E.;  Trushin, S. A. Chem. Phys. 1998, 232, 161‐174.  [20]. Robb, M. A.; Garavelli, M.; Olivucci, M.; Bernardi, F. Rev. Comput. Chem.  2000; Vol. 15, p 87‐146.  [21]. Toniolo,  A.;  Ben‐Nun, M.;  Martinez,  T.  J.  J. Phys.  Chem. A 2002,  106,  4679‐4689.  [22]. Yarkony, D. R. Faraday Discuss. 2004, 127, 325‐336.  [23]. Gilbert, A. T. B.; Besley, N. A.; Gill, P. M. W. J. Phys. Chem. A 2008, 112,  13164‐13171.  [24]. Frank,  I.; Hutter,  J.; Marx, D.; Parrinello, M.  J. Chem. Phys. 1998, 108,  4060‐4069.  [25]. Ziegler, T.; Rauk, A.; Baerends, E.  J. Theor. Chim. Acta 1977, 43, 261‐ 271.  [26]. Billeter, S. R.; Egli, D. J. Chem. Phys. 2006, 125, 224103(18).  [27]. Doltsinis, N. L.; Marx, D. Phys. Rev. Lett. 2002, 88, 166402‐4.  [28]. Runge, E.; Gross, E. K. U. Phys. Rev. Lett. 1984, 52, 997‐1000.  [29]. Levine, B. G.; Ko, C.; Quenneville, J.; Martinez, T. J. Mol. Phys. 2006, 104,  1039‐1051.  [30]. Mikhailov,  I. A.; Belfield, K. D.; Masunov, A. E.  J. Phys. Chem. A 2009,  113, 7080‐7089.  [31]. Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133‐A1138.  [32]. Adamo, C.; di Matteo, A.; Barone, V. Adv. Quantum Chem. 2000, 36, 45‐ 75.  [33]. Gritsenko, O. V.; Baerends, E. J. Theor. Chem. Acc. 1997, 96, 44‐50.  [34]. Takeda, R.;  Yamanaka,  S.;  Yamaguchi, K.  Int.  J. Quantum Chem 2005,  101, 658‐665.  [35]. Becke, A. D. J. Chem. Phys. 1997, 107, 8554‐8560.  [36]. Grimme, S.; Waletzke, M. J. Chem. Phys. 1999, 111, 5645‐5655.  [37]. Gimme, S. J. Chem. Phys. 2006, 124, 034108(16).  [38]. Yamanaka,  S.;  Nakata,  K.;  Takada,  T.;  Kusakabe,  K.;  Ugaide,  J.  M.;  Yamaguchi, K. Chem. Lett. 2006, 35, 242‐247.  [39]. Gutle, C.; Savin, A. Phys. Rev. A: At. Mol. Opt. Phys. 2007, 75, 032519‐17.  [40]. Langhoff,  P. W.;  Epstein,  S.  T.;  Karplus, M. Rev. Mod. Phys. 1972, 44,  602‐644.  [41]. Hansen, A. E.; Bouman, T. D. Mol. Phys. 1979, 37, 1713‐1724.  [42]. Casida, M. E. Recent Advances in Density­Functional Methods; Chong, D.  A., Ed. 1995; Vol. 3 of Part I.  [43]. Stratmann, R. E.; Scuseria, G. E.; Frisch, M. J. J. Chem. Phys. 1998, 109,  8218‐8224.  [44]. Dunning, T. H.; McKoy, V. J. Chem. Phys. 1967, 47, 1735‐1747.  [45]. Hirata, S.; Head‐Gordon, M. Chem. Phys. Lett. 1999, 314, 291‐299.  [46]. Savin, A.; Umrigar, C. J.; Gonze, X. Chem. Phys. Lett. 1998, 288, 391‐395.  160  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161    [47]. Neugebauer,  J.;  Baerends,  E.  J.;  Nooijen, M.  J. Chem. Phys. 2004, 121,  6155‐6166.  [48]. Hsu, C. P.; Hirata, S.; Head‐Gordon, M. J. Phys. Chem. A 2001, 105, 451‐ 458.  [49]. Hirata, S.; Head‐Gordon, M. Chem. Phys. Lett. 1999, 302, 375‐382.  [50]. Catalan, J.; de Paz, J. L. G. J. Chem. Phys. 2004, 120, 1864‐1872.  [51]. Starcke,  J.  H.; Wormit, M.;  Schirmer,  J.;  Dreuw, A. Chem. Phys. 2006,  329, 39‐49.  [52]. Wanko, M.; Garavelli, M.; Bernardi, F.; Niehaus, T. A.; Frauenheim, T.;  Elstner, M. J. Chem. Phys. 2004, 120, 1674‐1692.  [53]. Fantacci, S.; Migani, A.; Olivucci, M.  J. Phys. Chem. A 2004, 108, 1208‐ 1213.  [54]. Shao,  Y.  H.;  Head‐Gordon, M.;  Krylov,  A.  I.  J.  Chem. Phys. 2003,  118,  4807‐4818.  [55]. Guan,  J.  G.;  Wang,  F.;  Ziegler,  T.;  Cox,  H.  J.  Chem.  Phys.  2006,  125,  044314(9).  [56]. Maitra, N. T.; Zhang, F.; Cave, R. J.; Burke, K. J. Chem. Phys. 2004, 120,  5932‐5937.  [57]. Shibuya, T.; Rose, J.; McKoy, V. J. Chem. Phys. 1973, 58, 500‐507.  [58]. Jorgensen,  P.;  Swanstrom,  P.;  Yeager,  D.  L.;  Olsen,  J.  Int.  J. Quantum  Chem 1983, 23, 959‐971.  [59]. Hirata, S.; Nooijen, M.; Grabowski, I.; Bartlett, R. J. J. Chem. Phys. 2001,  114, 3919‐3928.  [60]. Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 7029‐7039.  [61]. DelBene, J. E.; Watts, J. D.; Bartlett, R. J. J. Chem. Phys. 1997, 106, 6051‐ 6060.  [62]. Dalgaard, E. J. Chem. Phys. 1980, 72, 816‐823.  [63]. Olsen, J.; Jorgensen, P. J. Chem. Phys. 1985, 82, 3235‐3264.  [64]. Vahtras, O.; Agren, H.; Jorgensen, P.; Jorgen, H.; Jensen, A.; Helgaker, T.;  Olsen, J. J. Chem. Phys. 1992, 97, 9178‐9187.  [65]. Sasagane, K.; Aiga, F.; Itoh, R. J. Chem. Phys. 1993, 99, 3738‐3778.  [66]. Luo, Y.; Agren, H.; Stafstrom, S. J. Phys. Chem. 1994, 98, 7782‐7789.  [67]. Salek, P.; Vahtras, O.; Guo, J. D.; Luo, Y.; Helgaker, T.; Agren, H. Chem.  Phys. Lett. 2003, 374, 446‐452.  [68]. Dalgaard, E.; Monkhorst, H. J. Phys. Rev. A: At. Mol. Opt. Phys. 1983, 28,  1217‐1222.  [69]. Larsen,  H.;  Jorgensen,  P.;  Olsen,  J.;  Helgaker,  T.  J.  Chem. Phys. 2000,  113, 8908‐8917.  [70]. Furche, F. J. Chem. Phys. 2001, 114, 5982‐5992.  [71]. Salek, P.; Vahtras, O.; Helgaker, T.; Agren, H. J. Chem. Phys. 2002, 117,  9630‐9645.  [72]. Cronstrand, P.; Luo, Y.; Agren, H. Adv. Quantum Chem. 2005, 50, 1‐21.  [73]. Knoester,  J.;  Mukamel,  S.  Phys.  Rev.  A:  At. Mol.  Opt.  Phys.  1989,  39,  1899‐1914.  [74]. Tretiak, S.; Mukamel, S. Chem. Rev. 2002, 102, 3171‐3212.  [75]. Tretiak, S.; Chernyak, V. J. Chem. Phys. 2003, 119, 8809‐8823.  [76]. Ziegler, T.; Seth, M.; Krykunov, M.; Autschbach, J.; Wang, F. The Journal  of Chemical Physics 2009, 130, 154102‐8.  [77]. Ziegler, T.; Seth, M.; Krykunov, M.; Autschbach, J.; Wang, F. Journal of  Molecular Structure: Theochem 2009, 914, 106‐109.  [78]. Katz, H. E. Chem. Mater. 2004, 16, 4748‐4756.  [79]. Dalton, L. R. J. Phys.: Condens. Matter 2003, 15, R897‐R934.  [80]. Perry, J. W.; Marder, S. R.; Meyers, F.; Lu, D.; Chen, G.; Goddard, W. A.;  Bredas, J. L.; Pierce, B. M. Polymers for Second­Order Nonlinear Optics,;  American Chemical Society: Washington, DC, 1995, 601, 45‐56.  [81]. Suponitsky, K. Y.; Timofeeva, T. V.; Antipin, M. Y. Usp. Khim. 2006, 75,  515‐556.  [82]. Schulten, K.; Karplus, M. Chem. Phys. Lett. 1972, 14, 305‐309.  [83]. Olchawa, R. Phys. B (Amsterdam, Neth.) 2000, 291, 29‐33.  [84]. Barford,  W.;  Bursill,  R.  J.;  Lavrentiev,  M.  Y.  Phys.  Rev.  B:  Condens.  Matter 2001, 63, 195108‐8.  [85]. Orr, B. J.; Ward, J. F. Mol. Phys. 1971, 20, 513‐526.  [86]. Wu, J. W.; Heflin, J. R.; Norwood, R. A.; Wong, K. Y.; Zamanikhamiri, O.;  Garito, A. F.; Kalyanaraman, P.; Sounik, J.  J. Opt. Soc. Am. B: Opt. Phys.  1989, 6, 707‐720.  [87]. McWilliams, P. C. M.; Hayden, G. W.; Soos, Z. G. Phys. Rev. B: Condens.  Matter 1991, 43, 9777‐9791.  [88]. Guo,  D.;  Mazumdar,  S.;  Dixit,  S.  N.;  Kajzar,  F.;  Jarka,  F.;  Kawabe,  Y.;  Peyghambarian, N. Phys. Rev. B: Condens. Matter 1993, 48, 1433‐1459.  [89]. Fitch, W. L.; McGregor, M.; Katritzky, A. R.; Lomaka, A.; Petrukhin, R.;  Karelson, M. J. Chem. Inf. Comput. Sci. 2002, 42, 830‐840.  [90]. Rumi, M.; Ehrlich, J. E.; Heikal, A. A.; Perry, J. W.; Barlow, S.; Hu, Z. Y.;  McCord‐Maughon,  D.;  Parker,  T.  C.;  Rockel,  H.;  Thayumanavan,  S.;  Marder,  S.  R.;  Beljonne, D.;  Bredas,  J.  L.  J. Am. Chem. Soc. 2000, 122,  9500‐9510.  [91]. Beljonne, D.; Shuai, Z.; SerranoAndres, L.; Bredas, J. L. Chem. Phys. Lett.  1997, 279, 1‐8.  [92]. Mikhailov, I. A.; Tafur, S.; Masunov, A. E. Phys. Rev. A: At. Mol. Opt. Phys.  2008, 77, 012510‐11.  [93]. Hashimoto, T.; Nakano, H.; Hirao, K.  J. Chem. Phys. 1996, 104,  6244‐ 6258.  [94]. Pariser, R. J. Chem. Phys. 1956, 24, 250‐268.  [95]. Hirao, K.; Nakano, H.; Nakayama, K.; Dupuis, M.  J. Chem. Phys. 1996,  105, 9227‐9239.  [96]. Malrieu, J. P.; Nebotgil, I.; Sanchezmarin, J. Pure Appl. Chem. 1984, 56,  1241‐1254.  [97]. Glushkov, V. N. Opt. Spectrosc. 2005, 99, 684‐689.  [98]. Löwdin, P.‐O.; Shull, H. Phys. Rev. 1956, 101, 1730.  [99]. Headgordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. Chem. Phys. Lett. 1994,  219, 21‐29.  [100]. Kurashige, Y.; Nakano, H.; Nakao, Y.; Hirao, K. Chem. Phys. Lett. 2004,  400, 425‐429.  [101]. Sekino, H.; Bartlett, R. J. Adv. Quantum Chem. 1999, 35, 149‐173.  [102]. Nakatsuji, H. Chem. Phys. Lett. 1991, 177, 331‐337.  [103]. Saha, B.; Ehara, M.; Nakatsuji, H. J. Chem. Phys. 2006, 125, 014316‐14.  [104]. Kitao, O.; Nakatsuji, H. Chem. Phys. Lett. 1988, 143, 528‐534.  [105]. Ostojic, B.; Domcke, W. Chem. Phys. 2001, 269, 1‐10.  [106]. Kawski, A.; Kuklinski, B.; Bojarski, P. Chem. Phys. 2006, 330, 307‐312.  [107]. Bublitz, G. U.; Boxer, S. G. Annu. Rev. Phys. Chem. 1997, 48, 213‐242.  [108]. Korter,  T. M.;  Borst,  D.  R.; Butler,  C.  J.;  Pratt,  D. W.  J. Am. Chem. Soc.  2001, 123, 96‐99.  [109]. Nakayama,  K.;  Nakano,  H.;  Hirao,  K.  Int.  J. Quantum Chem 1998, 66,  157‐175.  [110]. Dunning, T. H. J. Chem. Phys. 1989, 90, 1007‐1023.  [111]. Stanton,  J.  F.;  Gauss,  J.; Watts,  J.  D.;  Lauderdale, W.  J.;  Bartlett,  R.  J.;  Quantum  Theory  Project;  Departments  of  Chemistry  and  Physics,  University of Florida: Gainesville FL, 1993.  [112]. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.;  Cheeseman, J. R.; Montgomery, J., J. A. ; Vreven, T.; Kudin, K. N.; Burant,  J. C.; Millam,  J. M.;  Iyengar, S.  S.; Tomasi,  J.; Barone, V.; Mennucci, B.;  Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada,  M.;  Ehara,  M.;  Toyota,  K.;  Fukuda,  R.;  Hasegawa,  J.;  Ishida,  M.;  Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.;  Hratchian,  H.  P.;  Cross,  J.  B.;  Bakken,  V.;  Adamo,  C.;  Jaramillo,  J.;  Gomperts,  R.;  Stratmann,  R.  E.;  Yazyev,  O.;  Austin,  A.  J.;  Cammi,  R.;  Pomelli,  C.;  Ochterski,  J.  W.;  Ayala,  P.  Y.;  Morokuma,  K.;  Voth,  G.  A.;  Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels,  A.  D.;  Strain,  M.  C.;  Farkas,  O.;  Malick,  D.  K.;  Rabuck,  A.  D.;  Raghavachari,  K.;  Foresman,  J.  B.;  Ortiz,  J.  V.;  Cui,  Q.;  Baboul,  A.  G.;  Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz,  P.;  Komaromi,  I.;  Martin,  R.  L.;  Fox,  D.  J.;  Keith,  T.;  Al‐Laham, M.  A.;  Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.;  Johnson,  B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople,  J. A.; Revision E.1 ed.;  Gaussian, Inc.: Wallingford CT, 2004.  [113]. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.;  Cheeseman, J. R.; Montgomery, J., J. A. ; Vreven, T.; Kudin, K. N.; Burant,  J. C.; Millam,  J. M.;  Iyengar, S.  S.; Tomasi,  J.; Barone, V.; Mennucci, B.;  Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada,  M.;  Ehara,  M.;  Toyota,  K.;  Fukuda,  R.;  Hasegawa,  J.;  Ishida,  M.;  Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.;  Hratchian,  H.  P.;  Cross,  J.  B.;  Bakken,  V.;  Adamo,  C.;  Jaramillo,  J.;  Gomperts,  R.;  Stratmann,  R.  E.;  Yazyev,  O.;  Austin,  A.  J.;  Cammi,  R.;  Pomelli,  C.;  Ochterski,  J.  W.;  Ayala,  P.  Y.;  Morokuma,  K.;  Voth,  G.  A.;  Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels,  A.  D.;  Strain,  M.  C.;  Farkas,  O.;  Malick,  D.  K.;  Rabuck,  A.  D.;  Raghavachari,  K.;  Foresman,  J.  B.;  Ortiz,  J.  V.;  Cui,  Q.;  Baboul,  A.  G.;  Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz,  P.;  Komaromi,  I.;  Martin,  R.  L.;  Fox,  D.  J.;  Keith,  T.;  Al‐Laham, M.  A.;  Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.;  Johnson,  B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A.; Revision A.11 ed.;  Gaussian, Inc.: Pittsburgh PA, 1998.  [114]. Cumpston, B. H.; Ananthavel, S. P.; Barlow, S.; Dyer, D. L.; Ehrlich, J. E.;  Erskine,  L.  L.;  Heikal,  A.  A.;  Kuebler,  S.  M.;  Lee,  I.  Y.  S.;  McCord‐ Maughon, D.; Qin,  J. Q.; Rockel, H.; Rumi, M.; Wu, X. L.; Marder, S. R.;  Perry, J. W. Nature 1999, 398, 51‐54.  [115]. Kagotani, Y.; Miyajima, K.; Oohata, G.; Saito, S.; Ashida, M.; Edamatsu,  K.; Itoh, T. J. Lumin. 2005, 112, 113‐116.  [116]. Zipfel, W. R.; Williams, R. M.; Christie, R.; Nikitin, A. Y.; Hyman, B. T.;  Webb, W. W.  In National Academy of Sciences of  the United States of  America 2003, 100, 7075‐7080.  [117]. Wang,  C.  K.;  Macak,  P.;  Luo,  Y.;  Agren,  H.  J.  Chem.  Phys.  2001,  114,  9813‐9820.  [118]. Albota, M.; Beljonne, D.; Bredas, J. L.; Ehrlich, J. E.; Fu, J. Y.; Heikal, A. A.;  Hess, S. E.; Kogej, T.; Levin, M. D.; Marder, S. R.; McCord‐Maughon, D.;  Perry, J. W.; Rockel, H.; Rumi, M.; Subramaniam, C.; Webb, W. W.; Wu,  X. L.; Xu, C. Science 1998, 281, 1653‐1656.  [119]. Belfield, K. D.; Schafer, K. J.; Liu, Y. U.; Liu, J.; Ren, X. B.; Van Stryland, E.  W. J. Phys. Org. Chem. 2000, 13, 837‐849.  [120]. Belfield, K. D.; Morales, A. R.; Kang, B. S.; Hales, J. M.; Hagan, D. J.; Van  Stryland,  E.  W.;  Chapela,  V.  M.;  Percino,  J.  Chem.  Mater.  2004,  16,  4634‐4641.  [121]. Kim,  S.;  Wang,  Z.;  Hagan,  D.  J.;  Van  Stryland,  E.  W.;  Kobyakov,  A.;  Lederer, F.; Assanto, G. IEEE J. Quantum Electron. 1998, 34, 666‐672.  [122].  Ventelon,  L.;  Moreaux,  L.;  Mertz,  J.;  Blanchard‐Desce,  M.  Chem.  Commun. 1999, 2055‐2056.  [123]. Kim, O. K.; Lee, K. S.; Woo, H. Y.; Kim, K. S.; He, G. S.; Swiatkiewicz,  J.;  Prasad, P. N. Chem. Mater. 2000, 12, 284‐+.  [124]. Mikhailov,  I.  A.;  Bondar, M.  V.; Belfield,  K. D.; Masunov, A.  E.  J. Phys.  Chem. C 2009, 113, 20719‐20724.  Masunov and Mikhailov / European Journal of Chemistry 1 (2) (2010) 142­161  161  [125]. Suponitsky,  K.  Y.; Masunov,  A.  E.;  Antipin, M.  Y. Mendeleev Commun.  2009, 19, 311‐313.  [126]. Iordanov, T. D.; Davis, J. L.; Masunov, A. E.; Levenson, A.; Przhonska, O.  V.; Kachkovski, A. D. Int. J. Quantum Chem 2009, 109, 3592‐3601.  [127]. Suponitsky, K. Y.; Liao, Y.; Masunov, A. E.  J. Phys. Chem. A 2009, 113,  10994‐11001.  [128]. Toro, C.; De Boni, L.; Yao, S.; Ritchie, J. P.; Masunov, A. E.; Belfield, K. D.;  Hernandez, F. E. J. Chem. Phys. 2009, 130, 6.  [129]. Belfield,  K.  D.;  Bondar,  M.  V.;  Hernandez,  F.  E.;  Masunov,  A.  E.;  Mikhailov, I. A.; Morales, A. R.; Przhonska, O. V.; Yao, S. J. Phys. Chem. C  2009, 113, 4706‐4711.  [130]. Suponitsky,  K.  Y.; Masunov,  A.  E.;  Antipin, M.  Y. Mendeleev Commun.  2008, 18, 265‐267.  [131]. Suponitsky, K. Y.; Tafur, S.; Masunov, A. E. J. Chem. Phys. 2008, 129, 11.  [132]. Toro,  C.;  Thibert,  A.;  De  Boni,  L.; Masunov,  A.  E.;  Hernandez,  F.  E.  J.  Phys. Chem. B 2008, 112, 929‐937.  [133]. Kauffman, J. F.; Turner, J. M.; Alabugin, I. V.; Breiner, B.; Kovalenko, S.  V.; Badaeva, E. A.; Masunov, A.; Tretiak, S. J. Phys. Chem. A 2006, 110,  241‐251.  [134]. Masunov, A.; Tretiak, S.; Hong, J. W.; Liu, B.; Bazan, G. C. J. Chem. Phys.  2005, 122, 10.  [135]. Kobko, N.; Masunov, A.; Tretiak,  S. Chem. Phys. Lett. 2004, 392, 444‐ 451.  [136]. Fabian,  J.; Diaz, L. A.; Seifert, G.; Niehaus, T.  J. Mol. Struct. 2002, 594,  41‐53.  [137]. Hales, J. M.; Hagan, D. J.; Van Stryland, E. W.; Schafer, K. J.; Morales, A.  R.; Belfield, K. D.; Pacher, P.; Kwon, O.; Zojer, E.; Bredas,  J. L.  J. Chem.  Phys. 2004, 121, 3152‐3160.  [138]. Tafur, S.; Mikhailov, I.; Belfield, K.; Masunov, A. Computational Science  – LNCS 2009, 5545, 179‐188.  [139]. Belfield, K. D.; Bondar, M. V.; Hales, J. M.; Morales, A. R.; Przhonska, O.  V.; Schafer, K. J. J. Fluoresc. 2005, 15, 3‐11.  [140]. Belfield,  K.  D.;  Yao,  S.;  Morales,  A.  R.;  Hales,  J.  M.;  Hagan,  D.  J.;  Van  Stryland, E. W.; Chapela, V. M.; Percino,  J. Polym. Adv. Technol. 2005,  16, 150‐155.  [141]. Belfield,  K.  D.;  Bondar,  M.  V.;  Cohanoschi,  I.;  Hernandez,  F.  E.;  Kachkovsky, O. D.; Przhonska, O. V.; Yao, S. Appl. Opt. 2005, 44, 7232‐ 7238.  [142]. Schafer‐Hales, K.  J.; Belfield, K. D.; Yao, S.; Frederiksen, P. K.; Hales,  J.  M.; Kolattukudy, P. E. J. Biomed. Opt. 2005, 10.  [143]. Belfield, K. D.; Schafer, K. J. Chem. Mater. 2002, 14, 3656‐3662.  [144]. Day, P. N.; Nguyen, K. A.; Pachter, R. J. Phys. Chem. B 2005, 109, 1803‐ 1814.  [145]. Day,  P.  N.;  Nguyen,  K.  A.;  Pachter,  R.  J.  Chem.  Phys.  2006,  125,  094103(13).  [146]. Ohta, K.; Kamada, K. J. Chem. Phys. 2006, 124, 124303‐11.  [147]. Cronstrand, P.; Jansik, B.; Jonsson, D.; Luo, Y.; Agren, H. J. Chem. Phys.  2004, 121, 9239‐9246.  [148]. Masunov, A. M.; Tretiak, S. J. Phys. Chem. B 2004, 108, 899‐907.  [149]. Martin, R. L. J. Chem. Phys. 2003, 118, 4775‐4777.  [150]. Kokalj, A. Comput. Mater. Sci. 2003, 28, 155‐168. Code available from  http://www.xcrysden.org/   [151]. Belfield, K. D.; Bondar, M. V.; Hernandezt, F. E.; Przhonska, O. V.; Yao, S.  J. Phys. Chem. B 2007, 111, 12723‐12729.  [152]. Corredor, C. C. Ph. D. dissertation, University of Central Florida, 2007.  [153]. Gunnarsson, O.; Lundqvist, B. I. Phys. Rev. B: Condens. Matter 1976, 13,  4274‐4298.  [154]. Savin, A.; Colonna, F.;  Pollet, R.  Int.  J. Quantum Chem 2003, 93,  166‐ 190.  [155]. Della Sala, F.; Gorling, A. J. Chem. Phys. 2003, 118, 10439‐10454.  [156]. Moreira, I. D. R.; Costa, R.; Filatov, M.; Illas, F. J. Chem. Theory Comput.  2007, 3, 764‐774.  [157]. Cai, Z. L.; Reimers, J. R. J. Chem. Phys. 2000, 112, 527‐530.  [158]. Casida, M. E.; Ipatov, A. Abstr Pap Am Chem S 2006, 231, 1.  [159]. Cave, R.  J.;  Zhang, F.; Maitra, N. T.; Burke, K. Chem. Phys. Lett. 2004,  389, 39‐42.  [160]. Mikhailov,  I.  A.; Masunov,  A.  E. Computational  Science  ­ LNCS 2009,  5545, 169‐178.  [161]. Slater, J. C. Adv. Quantum Chem. 1972; Vol. 6, p 1‐92.  [162]. Noodleman, L.; Baerends, E. J. J. Am. Chem. Soc. 1984, 106, 2316‐2327.  [163]. Liberman, D. A. Phys. Rev. B: Condens. Matter 2000, 62, 6851‐6853.  [164]. Hu, C. P.; Sugino, O. J. Chem. Phys. 2007, 126, 074112(10).  [165]. Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, 2,  364‐382.  [166]. Takahashi, O.; Sumita, M. J. Mol. Struct. 2005, 731, 173‐175.  [167]. Schaftenaar, G.; Noordik, J. H. J. Comput. Aided Mol. Des. 2000, 14, 123‐ 134.  [168]. Kuthirummal,  N.;  Rudakov,  F. M.;  Evans,  C.  L.; Weber,  P. M.  J. Chem.  Phys. 2006, 125, 133307‐8.  [169]. Lochbrunner, S.; Fuss, W.; Schmid, W. E.; Kompa, K. L. J. Phys. Chem. A  1998, 102, 9334‐9344.  [170]. Woodward, R. B.; Hoffmann, R. J. Am. Chem. Soc. 1965, 87, 395‐397.  [171]. Longuet‐Higgins, H. C.; Abrahamson, E. W. J. Am. Chem. Soc. 1965, 87,  2045‐2046.  [172]. Michl, J. Molecular Photochemistry 1972, 4, 287‐314.  [173]. Celani, P.; Ottani, S.; Olivucci, M.; Bernardi, F.; Robb, M. A. J. Am. Chem.  Soc. 1994, 116, 10141‐10151.  [174]. Celani, P.; Bernardi,  F.; Robb, M. A.; Olivucci, M.  J. Phys. Chem. 1996,  100, 19364‐19366.  [175]. Garavelli, M.; Celani, P.; Fato, M.; Bearpark, M. J.; Smith, B. R.; Olivucci,  M.; Robb, M. A. J. Phys. Chem. A 1997, 101, 2023‐2032.  [176]. Tamura, H.; Nanbu, S.; Nakamura, H.; Ishida, T. Chem. Phys. Lett. 2005,  401, 487‐491.  [177]. Tamura, H.; Nanbu, S.; Ishida, T.; Nakamura, H. The Journal of Chemical  Physics 2006, 124, 084313‐13.  [178]. Guillaumont,  D.;  Kobayashi,  T.;  Kanda,  K.;  Miyasaka,  H.;  Uchida,  K.;  Kobatake, S.; Shibata, K.; Nakamura, S.; Irie, M. J. Phys. Chem. A 2002,  106, 7222‐7227.  [179]. Lafond,  C.;  Lessarda,  R.  A.;  Bolteb, M.;  Petkov,  I.  In Part  of  the  SPIE  Conference  on  Photopolymer  Device  Physics,  Chemistry,  and  Applications  IV.;  SPIE‐The  International  Society  for  Optical  Engineering, Vol.  3417: Québec, Canada, 1998, p 216‐227.  [180]. Andrasik, S. J.; Belfield, K. D.; Bondar, M. V.; Hernandez, F. E.; Morales,  A. R.; Przhonska, O. V.; Yao, S. ChemPhysChem 2007, 8, 399‐404.  [181]. Adronov,  A.;  Frechet,  J.  M.  J.;  He,  G.  S.;  Kim,  K.  S.;  Chung,  S.  J.;  Swiatkiewicz, J.; Prasad, P. N. Chem. Mater. 2000, 12, 2838‐2841.