Computer Simulation of the Nonlinear Optical Properties of Langmuir

Article pubs.acs.org/JPCC

Computer Simulation of the Nonlinear Optical Properties of Langmuir−Blodgett Films of a Squaraine Derivative G. Megariotis, A. Avramopoulos, M. G. Papadopoulos, and H. Reis* Institute of Biology, Medicinal Chemistry and Biotechnology, National Hellenic Research Foundation, Vassileos Konstandinou Avenue 48, 11635 Athens, Greece ABSTRACT: The linear and nonlinear optical properties (L&NLO) of a monolayer of (4-(N-methyl-N-(carboxypropyl)amino)-phenyl-4′-(N,Ndibutylamino)phenylsquaraine (SBA) were computed by a combination of molecular dynamics simulations, molecular quantum-mechanical, and discrete local-field model calculations to compute the macroscopic optical susceptibilities. Langmuir films on water and graphite were simulated and several structural properties investigated. In addition, SBA in water and dimethylsulfoxide (DMSO) solutions were simulated, and the aggregation behavior of SBA in these solvents was analyzed. The molecular (hyper)polarizabilities of SBA were calculated with density functional theory and validated with MP2 correlated ab initio calculations. An unusual large component perpendicular to the molecular dipole moment is obtained for the molecular first hyperpolarizability, which is a main contributor to the large in-plane anisotropy of the computed macroscopic second-harmonic generation signal of the Langmuir−Blodgett film. In order to circumvent the polarization catastrophe in the discrete local field model, a Thole-type interaction model was applied. The excitation spectrum of SBA was computed both in the gas phase and in water solution using the polarizable continuum model (PCM). The absorption spectrum of the SBA monomer is reasonably well described by the PCM model, but it fails for the spectrum of the dimer. polarizabilities. We also extended a discrete local field model8,9 in order to be able to predict the L&NLO susceptibilities of Langmuir−Blodgett films and applied this model to the simulated LB films. We chose to use 4-(N-methyl-N(carboxypropyl)amino)-phenyl-4′-(N,N-dibutylamino)phenylsquaraine (SBA, Figure 1) as our model system, which

1. INTRODUCTION The electronic structure and intermolecular interactions of squaraine dyes have been an area of intense interest for more than three decades.1−7 Squaraines can be described as quadrupolar molecules, which are highly polarizable and thus may exhibit large and unusual linear and nonlinear optical (L&NLO) properties. Apart from interesting molecular properties, such as negative second hyperpolarizabilities1,2 and negative electric-field-induced second-harmonic generation signals,3,4 there are strong intermolecular interactions between squaraine molecules in larger aggregates which can lead to quite unusal L&NLO properties. For example, Ashwell et al. observed strong second-harmonic generation (SHG) signals from Langmuir−Blodgett films composed of essentially centrosymmetric squaraine dyes.5 Molecular centrosymmetry excludes the possibility of dipolar SHG, and the observed signal was explained by formation of noncentrosymmetric aggregates, together with the occurrence of intermolecular charge transfer contributing to the SHG signal of the films. A fully computational treatment of the formation and properties of squaraine aggregates in matter would be highly desirable in order to investigate the unusual properties in more detail. Such a treatment will have to overcome several challenges in order to be successful: a reliable molecular description of the molecular properties has to be combined with an equally reliable description of the intermolecular interactions in the films. Here we report on some first steps in this direction. We used classical molecular dynamics to simulate thin films of a squaraine dye and performed quantum-chemical calculations for characterization of the molecular (hyper)© 2012 American Chemical Society

Figure 1. Schematic representation of a SBA molecule.

has been investigated extensively both in solution and in Langmuir and Langmuir−Blodgett films by Chen et al.6,7 Their results can thus be used as a valuable reference for theoretical predictions, in particular, for simulations of the macroscopic assemblies. Chen et al. found that the absorption spectrum of SBA in water is concentration dependent. This could quantitatively be interpreted by an equilibrium reaction between a dimeric form of SBA and the monomer. A H-dimer or cardpack association was assumed, concurrent with a considerable hypsochromic shift of the dimeric form. To test if our force fields are able to predict this kind of squaraine aggregation, we performed molecular simulations on solutions of SBA in water and DMSO Received: May 2, 2012 Revised: June 26, 2012 Published: July 17, 2012 15449

dx.doi.org/10.1021/jp304246a | J. Phys. Chem. C 2012, 116, 15449−15457

The Journal of Physical Chemistry C

Article

tions were done both in the gas phase and in the presence of solvents, applying the polarizable continuum model (PCM).38 In order to take into account the different conditions for electronic and nuclear solvation as accurately as possible, the state-specific solvation model of the vertical excitation as described by Improta et al.39 and implemented in Gaussian0933 was applied. All reported DFT computations were done at the B3LYP/6-31+G(d,p)-optimized geometry. 2.3. Local Field. In the discrete local field approximation a molecule is affected by local fields due to permanent and induced multipole moments of the molecules of the surrounding environment. Restricting treatment to the first nonvanishing multipole moment, the dipole, the permanent local field F̲0k at molecule k due the surrounding molecules k′ is

as well as ab initio calculations of the excitation spectra of the monomer and a dimer form of SBA both in the gas phase as well as in solution using the polarizable continuum model (PCM).

2. METHODS 2.1. Molecular Dynamics Simulations. All simulations were performed with the MD package GROMACS 4.0.3.10−12 Equations of motion were integrated with a time step of 1 fs. The NVT ensemble was employed for simulations of monolayers, and the temperature was kept constant at 300 K using the Nose−Hoover thermostat13,14 with a time constant of 1.0 ps. Solutions of SBA in water were simulated in the NPT ensemble (isotropic coupling) using the Parrinello−Rahman barostat15,16 with a time constant of 1.0 ps. In the case of the Langmuir film, the simulation box consists of two interfaces with a layer of water in between. Standard periodic boundary conditions were applied in all directions; effects of periodicity in the c direction (normal to the monolayers) were minimized using a large length of the box along this direction (∼3.2 nm, measured from the outermost carbon atoms of the alkyl chains). The dimensions of the box were a = 5.6734 nm, b = 6.79175 nm, c = 14.72388 nm. The approach adopted here for simulation of monolayers has been used generally in other studies of such systems.17−22 In the starting configuration, SBA molecules were placed at the two interfaces with the polar groups toward the water phase. The molecules were distributed randomly in the a−b plane, with the polar groups immersed in the water phase. After building the initial configuration of the system, a steepest descent energy minimization procedure was employed, followed by an equilibration MD run of 10 ns. Finally, a MD production run of 10 ns produced the trajectories which were then used for analysis. All these MD runs simulated a NVT ensemble. The topology file for the SBA was prepared with the PRODGR server23 using the GROMOS force field.24 All bonds were constrained to their equilibrium length with the LINCS algorithm,25 while the Ryckaert−Bellemans potential26,27 was applied for the dihedral angles of the SBA alkyl tails. The partial charges of SBA were calculated by fitting the electrostatic potential of a RHF/6-31G** wave function of the B3LYP/631G**-optimized SBA molecule using the CHelpG method introduced by Breneman and Wiberg.28 During the simulations, electrostatic interactions were treated with the Particle Mesh Ewald method.12,29,30 Water molecules were simulated by the SPC31 model in the case of the Langmuir film and SPC, SPC/ E32 models for solutions. 2.2. Molecular L&NLO Properties. The molecular static (hyper)polarizabilities of SBA were calculated with density functional theory (DFT) as well as correlated second-order Møller−Plesset theory (MP2) using GAUSSIAN 09.33 The B3LYP functional was used in the DFT computations using the aug-cc-pVDZ34 and 6-31+G(d,p)35 basis sets. The DFT (hyper)polarizabilities were computed either analytically or numerically by the finite field method. In the case of finite field calculations, the Romberg approach was applied in order to remove the higher order contaminations36 using a number of field strengths of magnitude 2mF, where m = 1, 2, 3, 4 and F = 0.0005 au. The excitation spectrum of the monomer and a dimer of SBA molecule was calculated with TDDFT, employing the B3LYP functional and the 6-31+G(d,p) basis set. TDDFT is often able to semiquantitatively describe lowlying electronic transitions in organic molecules.37 Computa-

N

F̲ ka0 =

∑

T (2)kk′, αβ( r̲ kk ′) · μ̲ keff ′β

k ′≠ k

where N is the number of molecules in the simulation box and Einstein summation is assumed for Cartesian indices (α,β). T(2) is the dipole field tensor. The effective dipole is given by μ̲ keff = μ̲ k0′ α + α k0′ , αβ · F̲ k0′ , α ′α

where α 0k′,αβ are components of the static f ree molecule dipole− dipole polarizability tensor. A more detailed description on the induction model is given in refs40 and 41. 2.4. Macroscopic Susceptibilities. The linear ( χ (1)(ω)) and second-harmonic generation (SHG, χ (2)(−2ω;ω,ω)) susceptibilities are computed according to 1 ε0v

χ (1) (ω) =

∑ α k(ω)· d k(ω) k

(1)

χ (2) ( −2ω; ω , ω) =

1 2ε0v

∑ d kT(2ω)· β k

k

( −2ω: ω , ω): d k(ω) d k(ω) (2)

where ⟨...⟩ indicates time averaging over the simulation trajectory, k labels the N molecules in the molecular simulation box with volume v, d k is a local field factor tensor,8 α k(ω), β k(−2ω: ω,ω) are the effective polarizability and first hyperpolarizability of molecule k in the phase, and ε0 is the permittivity of vacuum. The point-molecule multipole approximation is inadequate for the electrostatic response of molecules as large as SBA in condensed phases. For a semiquantitative description, we resort to the so-called submolecule treatment,42 which divides the molecule into a set of equivalent point submolecules over which the electric properties are equally distributed. Despite its simplicity, previous applications using this treatment to a variety of different molecular systems9,40,41,43,43 have shown that it can be surprisingly effective. After several tests with different partition models, we decided to use a model where a submolecule was assigned to each non-hydrogen atom, leading to 35 submolecules per molecule SBA. Application of the submolecule model is simple: the molecular properties α and β in eqs 1 and 2 are divided by the number of submolecules, and 15450



Article

χii(1)(ω = 0.0478 au), obtained with the second set α(ω = 0.0478 au), differed by a large factor from χii(1)(ω = 0.0239 au) (we choose, a bit arbitrarily, the value 2 as the limit), a polarization catastrophe was likely the reason. Thus, the parameter a was determined as the lowest value where for all snapshots χ(1)ii(ω = 0.0239 au)/2 < χ(1)ii(ω = 0.0478 au) < χ(1)ii(ω = 0.0239 au) × 2. This automatically includes the case of negative diagonal components. Clearly, it was a lucky coincidence that the two polarizabilities were just on the two opposite sides of the limit, where the catastrophe occurs. To treat the case where all polarizability tensors give rise to the catastrophe, one may first downscale the smallest polarizability tensors in such a way that no catastrophe occurs in any snapshot and use this downscaled tensor as a reference tensor for the method described above. The final value of the optimal parameter depends slightly on the density of the films, which differs considerably between the two simulations; defining a surface density σsurf by dividing the number of SBA molecules by the surface of the substrate, the value of a for the LB film on water with σsurf = 1.791 nm−2 is a = 1.1 and a = 0.8 for the Langmuir film on graphite, where σsurf = 3.326 nm−2. The values for a obtained here are considerably larger than those published in the literature based on the Thole model. This is probably due to the treatment of intramolecular interactions as of electrostatic origin in the Thole model,47 while in our model these interactions are included in the quantum-mechanical calculation of the reference structure. Intramolecular distances are generally smaller than intermolecular distances, and electrostatic interactions over the former require thus a stronger smearing.

the summations over the molecules are extended over all submolecules.40 In addition to the inadequacy of the one-point multipole approximation, the system investigated here suffers also from the so-called ‘polarization catastrophe’,45,46 as indicated by extremely large or even negative diagonal elements of the linear susceptibility tensor. In order to treat this problem, we chose to apply an adaptation of the “smearing” procedure first introduced by Thole in his dipole interaction model47 and later applied and refined, e.g., by van Duijnen and Swart48 and Ponder et al.49,50 Instead of calculating the interaction between point multipoles, this procedure computes the interaction between spatially distributed multipole densities, thus simulating in a more realistic manner the real electronic distribution in a molecule. The distribution function used to create smearedout densities should fulfill certain general conditions47 but is otherwise arbitrary. As shown by Burnham et al.,51 the electrostatic interaction between two smeared-out charge densities can be expressed by point multipoles interacting via a modified interaction tensor. This description is more convenient here, where the local field produced by the multipoles of the surrounding molecules is the relevant quantity. The local field effect is taken into account by the local field factor tensor d k in eqs 1 and 2, which is a 3N × 3 −1 , where matrix with components dki,j = ∑k′[ I − (εov)−1 L · α ]ki,k′j I is the 3N × 3N unit matrix, α is a 3N × 3N supermatrix with components αki,k′j = δkk′αij, and L is the Lorentz factor tensor, which is calculated by Ewald summation, using specifically the procedure outlined in ref 52. The modifications required to incorporate the smearing function were taken from ref 50. Following again ref 51 the smearing function, which is a shortrange correction, is only applied to the real space part of the Ewald sum, as the reciprocal sum accounts for long-range interactions. The functional form of the modified charge distribution used here is the same as the one applied previously, e.g., by Burnham et al.51 and Ponder et al.49,50 ρ( r ̲ ) =

⎛ | r̲ ij|3 ⎞ 3a ⎜ ⎟ exp −a 4π ⎜⎝ (αiαj)1/2 ⎟⎠

3. RESULTS AND DISCUSSION 3.1. Langmuir Films. The system under study contains 138 SBA and 4336 water molecules, providing 17 976 atoms in the simulation. The trajectory of 10 ns obtained from the MD simulations was used for analysis. First, the density profiles of water and SBA along the normal to interfaces were computed (Figure 2). The figure clearly shows the existence of two monolayers and the water phase.

(3)

where αi and αj are the mean polarizabilities of the submolecules at ri and rj, respectively, a is a parameter determining the width of the smeared charge distribution, and rij = ri − rj. In previous applications of Thole’s method, in which also all intramolecular interactions between atomic dipoles were modeled as electrostatic interactions, this parameter could be determined, together with the atomic polarizabilities, by fitting model molecular polarizabilities to either experimental or ab initio molecular polarizabilities, see, e.g., ref 48. This is not possible in the submolecule model, where intramolecular interactions are not modeled electrostatically, but taken as given by the computed molecular polarizability. Thus, to determine the parameter a, we used the fact that the occurrence of the polarization catastrophe is a consequence of the magnitude of the polarizabilities. It turned out that no polarization catastrophe occurred in any of the snapshots of the trajectories with the smaller set of the two frequencydependent polarizabilities, α(ω = 0.0239 au), leading to reasonable values for the instantaneous diagonal linear susceptibilities χii(1)(ω = 0.0239 au) for all snapshots. If for the same snapshot the corresponding instantaneous value of

Figure 2. Density profiles of water and SBA molecules along the normal to the monolayers.

As a first measure for the ordering of SBA, the Saupe ordering tensor Q 53 is computed by the following equation Q =

1 N

N

∑ i=1

{ 23 û û − 12 1} i i

(4)

where N is the number of SBA molecules, u î s the unit vector along the long axis of the rigid part of the ith molecule and 1 is 15451



Article

the unitary tensor of second order. The (dyadic) tensor Q is obtained as an instantaneous value for each configuration. As a criterion for the ordering we will consider the order parameter S, defined as the largest eigenvalue of Q . A more detailed analysis of this methodology can be found elsewhere.54−57 In our system, the SBA molecules will be oriented, on average, along a direction specified by a unit vector called the director, n. The director is defined here as that eigenvector of Q which corresponds to its largest eigenvalue. The order parameter S for our system is 0.518 ± 0.004. For an isotropic phase, values close to zero would be expected. Such a case is a common liquid with its molecules moving without preferential orientation. As regards the orientation of the director, this stays practically unchanged during the simulations, as indicated by the normalized autocorrelation function57 C(t ) = ⟨ n̲ (0) · n̲ (t )⟩

values are shown in Table 1 for the two alkyl tails of the SBA molecules. The values are close to zero, which indicates that Table 1. Deuterium Order Parameter SCD for the Alkyl Tail (Langmuir films)

Pxx + Pyy ⎫ Lz ⎧ ⎨Pzz − ⎬ 2 n⎩ ⎭

(5)

(6)

where Lz is the height of the box and n the number of surfaces (in our case n = 2). The Pii are components of the pressure tensor along the unit cell axis i. The value for γ obtained from our simulation is 62.5 ± 5.5 mN/m. This is in reasonable agreement with the experimental value of 67.0 mN/m, estimated from the experimental surface pressure π of ∼5 mN/m6 by γ = γwater − π and γwater = 72 mN/m at room temperature.59 Defining a molecular axis to be along the line of largest extension of the nearly rigid part of the molecule between the two nitrogen atoms, we computed the average tilt angle of this molecular axis of SBA with respect to the c axis. The obtained value was 51.9° ± 2.1°, which is in good agreement with the tilt angle of 50° ± 2° reported in ref 6, estimated by analyzing the difference between s- and p-linear polarized light of a monolayer of SBA on quartz. This comparison relies on the transition moment being along the molecular axis, as defined above. As shown below in section 3.5, this is in fact borne out by the ab initio calculations. A structural property that provides useful information about the structural order in the alkyl chains is connected with the molecular order tensor parameters Sij, defined by Sij =

1 ⟨3 cos θi cos θj − δij⟩ 2

alkyl tail 2

0.03 0.01

0.02 0.02

there is no preferential orientation. Unfortunately, no experimental results are available for comparison. The mean area per molecule (mA/M) in the simulation was 55.83 Å2, very close to the first maximum of the experimental surface pressure−area isotherm reported in ref 6. We also tried to simulate monolayers with different values for mA/M but were only partially successful: apart from simulations with mA/ M from 60 to 55 Å2, all films with different mA/M values collapsed. 3.2. Langmuir−Blodgett Film of SBA. A monolayer of SBA was also simulated on a solid substrate, specifically graphite. As this is a hydrophobic substrate, the alkyl tails of SBA are in the vicinity to the graphite atoms whereas the polar part is pointing away from the surface. The substrate is simulated by three graphene planes which are frozen during the simulation and interact with the SBA layer only through van der Waals interactions. Following a previous study of a similar simulation,61 the carbon atoms of graphite were simulated using the atoms of type CB of the GROMOS force field. The dimensions of the surface planes are 4.49 and 5.34 nm in the a and b direction, respectively, and 7.82 nm in the c direction. In this case the vacuum layer added in order to minimize interlayer interactions in the c direction is about 4 nm. The monolayer of SBA consists of 79 molecules, thus yielding a higher surface density than in the Langmuir film on water. In fact, the mA/M of this simulation is 30.1 Å2, thus smaller than the experimentally estimated limiting mA/M of 38 Å2 for SBA monolayers on water.6 Trials to simulate films with larger mA/ M values were unfortunately unsuccessful. The order parameter S of a 10 ns trajectory is equal to 0.600 ± 0.002, considerably larger than the corresponding value of the Langmuir film, probably as a consequence of the larger density of the film on graphite. The normalized autocorrelation function of the director (see eq 2) is very close to 1, never decreasing below 0.99. Again, the SCD values for the alkyl tails are close to zero, similar as in the case of the Langmuir film (see previous subsection and Table 1). Additionally, the tilt of the SBA molecular axis with respect to the c axis was found to be 41.3° ± 1.2°, slightly smaller than for the Langmuir film. 3.3. Solution of SBA in Water. The initial configurations for the simulations of solutions of SBA in water were constructed by placing SBA molecules randomly in a box which was then filled with water molecules. Four different concentrations of SBA were simulated, see Table 2. Two different force fields were used for water, i.e., SPC and SPC/E. The latter corrects for some of the deficiencies of the former, e.g., a greatly underestimated dielectric constant, see refs 60, 62, and 63. Table 2 shows the densities and dielectric constants of water−SBA solutions. The effect of increasing SBA concentration is a slight increase of the density and to decrease the dielectric constant. In order to investigate the possibility of SBA dimerization in water, we first computed the number of hydrogen bonds

In our simulations, the correlation function never falls below the value 0.94. Surface tension is an important experimental quantity for interfaces. When surfactants are added to the air−liquid interface, the surface tension is reduced due to formation of a monolayer. For molecular simulations as applied here, the surface tension (γ) can be computed from17,12,58 γ=

alkyl tail 1

(8)

where θi is the angle between the molecular axis i and the normal of the bilayer (c axis). The definition of the molecular axis for a CH2 group in an alkyl chain can be found in ref 60. Using values Saa and Sbb of a set of three united atoms of the chain, the quantity SCD = 2/3Saa + 1/3Sbb can be computed, which may be related to order parameters determined from 2H NMR experiments. Details about the computation of the properties Sij can be found elsewhere.12,61 For an alkyl chain of four carbons, two values for SCD can be defined; their computed 15452



Article

is 0.45 ± 0.02 nm and the average minimal distance between any atom of the two squarainyl cores is 0.32 ± 0.01 nm. In order to elucidate the structure of the dimer, we define two approximately perpendicular vectors in each SBA molecule, the first one pointing from the N atom of the N,N-dibutylamino group to the N atom of the N-methyl-N-carboxypropylamino group while the other vector points from one O atom to the other of the central squaryl unit. The computed average angle between the N−N vectors of the two molecules is 155° ± 10°, while the angle between the two O−O vectors is 30°± 5°. Thus, the two molecules are approximately oriented antiparallel to each other, with an angle of 30° between the long axes, and the planes of the molecules are approximately parallel to each other. As already mentioned, Chen et al.6 explained the concentration dependence of the absorption spectrum of SBA in water by a dimerization of SBA into H dimers. The simulated dimer structure found here in water with roughly antiparallel dipole moments is not the typical form of a H dimer, where the dipoles are generally aligned in parallel. Nevertheless, the strong ‘local’ C−O dipoles are still in an unfavorable parallel position, which then may lead to the blue shift observed for Hdimer structures. This question will be taken up again in section 3.5. Experimentally, the absorption spectrum of a solution of SBA in DMSO shows no concentration dependence and corresponds to the spectrum of the monomeric form.6 To test this we performed a simulation of SBA molecules in DMSO solvent with an average concentration of 0.0153 mol/L. In good accord with experiment, no dimerization was observed over the whole range of 20 ns production run (see Figure 3). We note that the lowest concentrations of the simulated solutions are much larger (∼0.02 mol/L) than those on which the absorption measurements were performed (∼10−6 mol/L6). 3.4. Molecular L&NLO Properties. In Table 4 the diagonal components of the static- and frequency-dependent

Table 2. Some Properties of Water−SBA Solutions water model

SPC

Nmol/unit cell

C (mol/ L)

ρ (kg/m3)

2 4 6 8

0.015 0.030 0.046 0.061

978.1 979.3 980.7 982.0

SPC/E ε

C (mol/ L)

ρ (kg/m3)

ε

63.5 63.8 64.5 64.5

0.015 0.031 0.047 0.062

999.6 1000.9 1002.2 1003.4

69.4 68.1 68.3 70.2

between SBA molecules. The standard Gromacs rules were applied: a hydrogen bond is assumed to be present if the distance between hydrogen and the acceptor atom is smaller than 0.35 nm and if the angle between acceptor−donor− hydrogen atoms is larger than 150°. As the values collected in Table 3 show, only the highest concentrated solution shows the occurrence of a sizable number of hydrogen bonds between SBA molecules. Table 3. Number NHB of Hydrogen Bonds Between SBA Molecules in the Simulated Solutions (normalized per SBA molecule) Nmol SBA/unit cell

SPC

SPC/E

2 4 6 8

0.000 0.005 0.015 0.300

0.000 0.152 0.227 0.334

However, visual inspection of the simulated structures in the solution with lowest concentration with SPC/E water shows that after about 8 ns a dimer of SBA molecules is formed and stays stable for the following 12 ns of the simulation (see Figure 3). Analyzing the last 10 ns of the simulation, the average distance of the center of the cyclobutadiene rings in the dimer

Figure 3. Temporal change of the distance r between the two closest SBA molecules in simulated DMSO and water solutions of lowest concentration, measured as the distance between the center of the cyclobutadiene rings, see text. 15453



Article

diagonal component βzzz, is reminiscent of the class of ‘twodimensional’ NLO chromophores,64 which are characterized by large nondiagonal components βzii, where z (≠i) is the dipole moment direction, and whose structure generally involves more than one charge-transfer pathways, arranged in a nonlinear way. Chen et al.3,4 previously reported that nonsymmetrically substituted squaraines in CHCl3 may show either a positive or a negative EFISH signal. In the case of SBA, the EFISH signal, which is essentially given by μz∑i=x,y,zβzii, is with 233 au positive (in the gasphase). That the dominant component of β is along a direction perpendicular to the dipole moment means that it will not contribute to the EFISH signal of a solution of SBA. On the other hand, it would contribute to the hyper-Rayleigh scattering (HRS) signal. Unfortunately, experimental EFISH or HRS values on SBA solutions, with which the ab initio predictions could be tested, have not been reported yet. Meyers et al.,65 using INDO calculations, have shown that the often invoked two-state approximation is not an adequate approximation for the first hyperpolarizability of several squaraine derivatives. This is also the case for SBA: Using the TDDFT/B3LYP-computed properties, the two-state contribution to the static βxxx, calculated according to

Table 4. Selected Components of the Dipole (μ), Electronic Polarizability (αii(ω)), and First (βjii(−2ω;ω,ω)) and Second (γiiii) Hyperpolarizabilities of SBA, Computed with the B3LYP Functional and Different Basis Sets (all values in au) aug-ccpVDZa

6-31+G** ω = 0 au μz αxx αyy αzz α̅ βxxx, βxyy, βxzz βyyy, βyxx, βyzz βzzz, βzxx, βzyy, βxyz γxxxx/103 γyyyy/103 γzzzz/103

−0.479 (−0.630b) 721.88 438.31 587.68 582.62 1368, −132, 255 34, −15, 205 48, −669, 135, −125 168 (170b) 122 (125b) 43 (46b)

ω = 0.0239 au

ω = 0 au −0.510 b

753.34 (753.56 ) 444.90 (444.86b) 609.46 (609.50b) 602.57 (602.64b) 1900, −199, 365 (1893, −206, 355)b 74, −106, 237 (90, −78, 252)b 17, −998, 205, −179 (5, −990, 206, −199)b

739.90 417.48 602.99 586.79 1345 43 −10 270 133 100

a

Only diagonal components were computed with the aug-cc-pVDZ basis. bEffective properties; local field (Fx, Fy, Fz) = (0.000,0.0002,0.0003) au.

βiiitwo‐state = 6

electronic polarizability and first as well as static second hyperpolarizabilities of the SBA molecule are shown. The two basis sets employed yield similar values for all components with the exception of the second hyperpolarizabilities, which are considerably larger for the aug-cc-pVDZ basis set, indicating that for these properties the smaller Pople basis set is not sufficient. To test the adequacy of the B3LYP functional for the electric properties of SBA, the diagonal components of α and β in x direction were also computed at the MP2/6-31+G(d,p) level of theory. The values obtained (αxx = 721.2 au, βxxx = 1700 au) are in reasonable agreement with those computed at the B3LYP/6-31+G(d,p) level (αxx = 721.9 au, βxxx = 1368 au). We mention that the second hyperpolarizabilities of SBA are not negative, in contrast to what has been found for several symmetrically substituted squaraine dyes.1,2 For the first hyperpolarizability, the diagonal element in the x direction, perpendicular to the dipole axis, is the dominant component (see Figure 4 for a visualization of the coordinate system) with the next largest component (βzxx) about a factor of 2 smaller. The comparatively large value for βzxx, together with the small

2 μ02, (μ − μ0, i ) i 2, i 2 ω02

where μ02,i is the ith component of the transition dipole from the ground state to the dominant low-lying dipole allowed excited state (which is the second excited state for SBA at the TDDFT/B3LYP level), μk,i is the ith component of the dipole in state k, and ω02 is the energy for the transition 0 → 2, is ∼−200 au, compared to 1368 au for the full value. 3.5. Electronic Excitations. For both the monomer molecule (Figure 4) and a dimer of SBA molecules with approximately parallel dipoles (Figure 5) the lowest lying dipole-allowed transition was calculated with TDDFT. A small blue shift (Table 5) of the first dipole-allowed excitation in going from the monomer (528.5) to the dimer (516.7) was found in the gas phase. The first dipole-allowed excitation band has also been computed in water (H2O) and dimethylsulfoxide (DMSO) solution using the polarizable continuum model

Figure 5. Structure of the optimized dimer; hydrogen atoms not shown.

Figure 4. Structure and coordinate system of the optimized monomer. 15454



Article

Table 5. Wavelength λ and Oscillator Strength f of the First Strongly Dipole-Allowed Excitation of the Monomer and Dimer SBA Computed in the Gas Phase and Solutions, Calculated with the B3LYP Functional and 6-31+G(d,p) Basis Set monomer λ (nm) exp.c f

dimer

gas opt.a

H2O opt.a

H2O simul.b,d

DMSO

gas

H2O opt.a

H2O simul.b

DMSO

528.5

639.4 650 2.21

558.9, 559.3

637.2 656 2.20

516.7

634.4 594 2.95

641, 642

632.1

1.59

2.08, 2.07

2.92

1.72, 2.70

3.04

a

Optimized monomer or optimized dimer structure (Figure 5). bStructure taken from the simulated trajectory of the lowest concentrated water solution and optimized with PCM/6-31+G**/B3LYP (Figure 6). cReference 6 dThe two entries refer to the two monomers of the dimer, optimized separately with PCM/6-31+G**/B3LYP.

(PCM). The solvent effect on the computed spectrum is remarkable, red shifting the gas phase values by about 120 nm, both for the monomer and for the dimer. The change of permanent dipole moment between the ground and the excited state cannot explain this strong solvatochromic effect: the dipole moment in the relevant excited state is with 0.490 au only marginally larger than the ground state dipole moment (0.479 au) at the B3LYP/6-31+G** level. Thus, either large changes of induced moments and/or higher multipoles are presumably responsible, which would be in accordance with the ‘quadrupolar’ character of the molecule. The computed values for monomer absorption in solution compare well with the corresponding experimental values in water and DMSO (see Table 5). The difference between the excitation in water and in DMSO is with 2.2 nm small, again in reasonable agreement with experiment,6 where a 6 nm difference was reported. However, the computed blue shift upon dimerization is reduced from the gas-phase value of 11.8 nm to ∼5 nm in both solvents. This is much smaller than the experimental value of 56 nm reported in ref 6, which was deduced from the concentration dependence of the absorption spectrum. As shown in section 3.3, the molecular simulation of a solution of SBA in water predicts a SBA dimer with quite a different geometry than the optimized geometry with parallel dipoles analyzed above. To check if this structural difference may have an effect on the calculated absorption spectrum, we extracted a typical dimer from the simulated trajectory, optimized it, as well as both monomers at the PCM/6-31+G**/B3LYP level, and finally computed the excitation parameters for the optimized dimer and the two monomers. Using the PCM model for water, the values shown in Table 5 were obtained. They show a considerable red shift upon dimerization. Two very close lying excitations are predicted for the dimer, both comparable to that of the optimized dimer, but the excitations of the two optimized monomers, taken from the dimer, are considerably higher than that of the optimized monomer. Overall, the values agree reasonably well with the experimental results but with the assignments of the bands to the dimer and monomer forms reversed. This is in contradiction to the concentration dependence of the absorption spectrum. We tentatively conclude from all the results that the PCM continuum model is not able to describe the effect of dimerization on the absorption spectrum of SBA in water accurately. Further investigations of this problem in the framework of the Förster exciton model66 are underway. We note finally that the transition moment of the monomer is directed along the long axis of the rigid part of the molecule, in agreement with experimental observations on a similar squarilium dye.67 3.6. Linear and SHG Susceptibilities of the Langmuir− Blodgett Films. Two model Langmuir−Blodgett films were

Figure 6. Structure of a dimer taken from the simulated trajectory of 2 SBA molecules in water and optimized at the PCM/6-31+G**/ B3LYP level; hydrogen atoms not shown for clarity.

created from the trajectories of the two simulations described in the previous sections, in order to compute their susceptibilities. The first approximation of a Langmuir−Blodgett film was produced from the results of the Langmuir-film simulation of SBA on a thin film of water by extracting the two monolayers on each side of the water−film separately and treating each as a model for a LB film of SBA. A second model of a LB film was created by removing the solid substrate layer in the simulation of SBA on graphite. To compute the effective molecular properties, which are required for calculation of the macroscopic susceptibilities, the local field F has to be known. It was computed according to the procedure outlined in section 2.3 for the LB film taken from the water simulation and found to be quite small, leading accordingly to small changes of the L&NLO properties (see Table 4, values in brackets). As usual,9,40,41,43 the effect of the local field is larger on the hyperpolarizabilities than on the polarizabilities. The effective properties obtained for this film were used for all subsequent susceptibility calculations. Table 6 shows the calculated linear and second-harmonic generation susceptibilities of the two films, computed according to eqs 3 and 4 and using the in-field, frequency-dependent properties given in Table 4. The coordinate system is oriented in such a way that the c axis is perpendicular to the surface of the original substrate (water or graphite). The numerical values computed for the susceptibilities itself are not meaningful, as the volume occurring in eqs 3 and 4 is that of the whole unit cell, thus taking into account the vacuum layer above the films, which needs to be thick enough to minimize ‘interlayer’ interactions, which would otherwise distort the single-layer susceptibilities we are interested in. For purposes of comparability, the thickness of the vacuum layer was kept equal (∼14.7 nm) for both films. 15455



Article

Table 6. Linear (χ(1)(ω)) and SHG (χ(2)(−2ω;ω,ω)/(pm V−1)) Susceptibilites of the Two Model LB Films Computed with the Smeared Density (eq 3) with a = 1.1 for the Water Film and a = 0.8 for the Film on Graphite; Statistical Errors are

Computer Simulation of the Nonlinear Optical Properties of Langmuir

Recommend Documents