4536
J. Phys. Chem. C 2010, 114, 4536–4543
Melting of an Anchored Bilayer: Molecular Dynamics Simulations of the Structural Transition in (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) Vikrant V. Naik and S. Vasudevan* Department of Inorganic and Physical Chemistry, Indian Institute of Science, Bangalore 560012, India ReceiVed: October 28, 2009; ReVised Manuscript ReceiVed: February 1, 2010
The thermally driven structural phase transition in the organic-inorganic hybrid perovskite (CnH2n+1NH3)2PbI4 has been investigated using molecular dynamics (MD) simulations. This system consists of positively charged alkyl-amine chains anchored to a rigid negatively charged PbI4 sheet with the chains organized as bilayers with a herringbone arrangement. Atomistic simulations were performed using an isothermal-isobaric ensemble over a wide temperature range from 65 to 665 K for different alkyl chain lengths, n ) 12, 14, 16, and 18. The simulations are able to reproduce the essential features of the experimental observations of this system, including the existence of a transition, the linear variation of the transition temperature with alkyl chain length, and the expansion of the bilayer thickness at the transition. By use of the distance fluctuation criteria, it is shown that the transition is associated with a melting of the alkyl chains of the anchored bilayer. An analysis of the conformation of the alkyl chains shows increased disorder in the form of gauche defects above the melting transition. Simulations also show that the melting transition is characterized by the complete disappearance of all-trans alkyl chains in the anchored bilayer, in agreement with experimental observations. A conformationally disordered chain has a larger effective cross-sectional area, and above the transition a uniformly tilted arrangement of the anchored chains can no longer be sustained. At the melt the angular distribution of the orientation of the chains are no longer uniform; the chains are splayed allowing for increased space for individual chains of the anchored bilayer. This is reflected in a sharp rise in the ratio of the mean head-to-head to tail-to-tail distance of the chains of the bilayer at the transition resulting in an expansion of the bilayer thickness. The present MD simulations provide a simple explanation as to how changes in conformation of individual alkyl-chains gives rise to the observed increase in the interlayer lattice spacing of (CnH2n+1NH3)2PbI4 at the melting transition. Introduction Bilayers formed by molecules that possess long alkyl hydrophobic tails are ubiquitous in the natural world manifesting both in biological systems as well as in chemistry. The lipid bilayer is an integral feature of cell membranes of living systems with functions that are of critical importance to the life of the cell. In chemistry they occur in vesicles,1 Langmuir-Blodgett films,2-4 supported bilayers,5,6 organoclays,7 self-assembled monolayers,3,8 surfactant intercalated solids,9 and in organicinorganic hybrid materials.10 Unlike free lipid bilayers, the bilayer in the latter examples, the surfactant intercalated solids, orgnaoclays, and organic-inorganic hybrids have one end anchored to an inorganic lattice that is quite rigid. These materials that possess an anchored or tethered bilayer have found utility as adsorbents of organic pollutants in soil and water remediation11-14 as well as in rheological control in paints and grease.15 They are also of interest in their own right as it has been shown both experimentally and from theoretical considerations that they can exhibit phase transitions and in analogy to similar transitions observed in lipid bilayers have been labeled as a melting transition.16-18 Many organic-inorganic hybrid systems based on the layered perovskite structure containing long (n > 12) alkylammonium chains are known to exhibit one or more phase transitions with temperature.19-25 Most studies have concluded that the main * To whom correspondence should be addressed. E-mail: svipc@ ipc.iisc.ernet.in. Phone: +91-80-2293-2661. Fax: +91-80-2360-1552/0683.
transition appearing at higher temperatures is a melting of the alkyl chains, a behavior similar to that in lipid bilayers.1 A notable difference is that unlike in lipid bilayers where individual molecules can undergo lateral diffusion and also flip-flop between layers the anchored bilayer is characterized by the total absence of translational mobility. The degrees of freedom of the alkyl chains of the anchored bilayer are restricted to changes in conformation. In lipid bilayers it is well established that the conformation and its dynamics plays a central role in defining the temperature dependence and the critical temperatures associated with the different phases of the bilayer. The anchored bilayer in these organic-inorganic hybrids is therefore of interest as a prototypical system for the melting behavior of the more complex lipid bilayers. They have the advantage that unlike lipid bilayers that are essentially fluidlike systems the organicinorganic hybrids are solids and permit investigations by the more robust solid-state spectroscopic techniques. Here we report investigations of a molecular dynamics (MD) simulations of the phase transition in the anchored bilayer of the organicinorganic hybrid (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) that consists of rigid PbI4 sheets flanked on either side by alkylammonium chains as depicted in Figure 1. This system exhibits a first-order endothermic transition accompanied by a dilation of the lattice and increased conformational disorder of the alkyl chains.19 A large number of MD simulation studies on the anchored bilayer in organoclays and surfactant intercalated layered solids have been reported, primarily focusing on the structure and
10.1021/jp910300v 2010 American Chemical Society Published on Web 02/18/2010
Melting of an Anchored Bilayer
J. Phys. Chem. C, Vol. 114, No. 10, 2010 4537
Figure 1. Snapshot of the simulation of (CnH2n+1NH3)2PbI4 (n ) 18) showing the alkylchains of the anchored bilayer in a herringbone arrangement. C is gray, H white, N blue, I magenta, and the PbI4 polyhedra green.
layering phenomena, either as function of the alkyl chain length or grafting density.26-31 MD simulations have been used to understand changes in conformation and basal plane spacings of alkyl chains anchored to mica sheets with temperature.26,32 At higher temperatures an increased basal spacing and gauche conformational disorder was observed. On the basis of these studies as well as experimental results, a relation between packing or anchoring densities and thermal transitions of the anchored bilayer was established.16 Atomistic MD simulations of the conformational mobility and dynamics of alkyl chains in the anchored bilayer of a surfactant intercalated solid have also been reported.33 These studies have shown how differences in the extent of conformational disorder as estimated by different spectroscopic techniques may be reconciled. A closely related system that has also been widely studied by MD techniques are the ultrathin films of alkanes physisorbed on graphite.34-37 Here too the focus has been on the conformation, structural arrangement, and the melting transition of the alkane chains. Unlike the earlier examples, however, lateral mobility is permitted, but as the chains are confined to the surface they do provide another example of the melting transition in twodimensions. The present MD simulations of the anchored bilayer in the PbI4 perovskite based organicsinorganic hybrid was prompted, in part, by a recent experimental report using vibrational spectroscopy that showed for all alkyl chain lengths the “melt” is characterized by the complete absence of all-trans chains.19,20 The objectives of the present study are 3-fold: First, to establish that atomistic MD simulations of this system do indeed exhibit a phase transition that satisfies or meets the criteria for being labeled a melting transition. Second, to investigate the thermal evolution of conformational disorder, establish the “molecularity” of the transition and see if the simulations can reproduce the experimental observation of the disappearance of planar alltrans chains at the transition. And finally, to establish how changes in the conformation of individual chains leads to the observed macroscopic changes in the interlayer lattice spacing of the anchored bilayer at the transition. Methodology Molecular dynamics simulations were performed using the Materials Studio suite38 running on an IBM M PRO Intellistation workstation. The first step in the simulation was the preparation of the unit cell of the (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) structure (henceforth abbreviated as Cn-PbI4). The unit cells were
constructed from that of (C4H9NH3)2PbI4 (C4-PbI4) for which the crystal data has been reported.39 This compound crystallizes in the Pbca space group with lattice parameters of a ) 8.8632 Å, b ) 8.6816 Å, c ) 27.570 Å, and R ) β ) γ ) 90°. The unit cells for (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) were constructed by first extending the c parameter of C4-PbI4 to the reported experimentally determined d value of the appropriate Cn-PbI4 compound while keeping the other lattice parameters a, b, R, β, and γ constant. In the second step the hydrocarbon chain of C4-PbI4 was extended by substituting the H atom of the terminal methyl group by an all trans hydrocarbon chain segment of appropriate length to obtain Cn-PbI4 with n ) 12, 14, 16, 18. Care was taken to ensure that all methyl groups in all chains were present in an all-trans conformation and that the chains were not interdigitated. The position of the headgroup was left unchanged from that of C4-PbI4 structure. The final lattice parameters of the unit cell were a ) 8.8632 Å, b ) 8.6816 Å, R ) β ) γ ) 90°, and the value of c the same as that obtained from X-ray diffraction. The unit cell thus obtained had the alkyl chains arranged in a herringbone pattern as shown in Figure 2 for C18-PbI4. A super cell was then constructed containing 64 such crystallographic unit cells. The super cell parameters were a ) 35.4528, b ) 34.7264, R ) β ) γ ) 90°, equivalent to 4 × 4 PbI4 unit cells in the ab plane and a two layer repeat in the c direction with an initial interlayer spacing of 33.6 Å for C18, 31.4 Å for C16, 28.2 Å for C14, and 26.1 Å for C12. The final superlattice formula was (CnH2n+1NH3)128Pb64I256. Except for the periodic boundary condition, no other symmetry constraints were imposed. The structure was treated as triclinic (P1), and all lattice parameters were treated as independent variables in the simulation. Charges were calculated separately for PbI4 layers and for the alkyl ammonium chains by assuming that the inorganic PbI4 layer of the unit cell carried a charge of -8e while each alkyl ammonium chain a unit positive charge. Charges on individual atoms of the PbI4 sheet were obtained by a DFT calculation using the DMol3 module of Material Studio.38 A double numerical with polarization (DNP) basis set with the PerdewWang (PW91) functional was used. Charges on atoms were obtained from a Mulliken population analysis using the eigenstates and eigenvectors from the above calculation. The charge on the Pb atom was 0.473. Iodine atoms at the corners of octahedra carried a negative charge of -0.778, while those sandwiched between two octahedra had a charge of -0.438. The alkylammonium chains have a unit positive charge. The
4538
J. Phys. Chem. C, Vol. 114, No. 10, 2010
Naik and Vasudevan region of the phase transition was identified, MD simulations were performed in steps of 10 K in this region. Results and Discussion
Figure 2. Unit cell of (C18H37NH3)2PbI4 as (a) viewed along the c axis and (b) viewed along the b axis.
charges on individual atoms of the chains were obtained by a Mulliken population analysis of the molecular orbitals calculated by the Hartree-Fock method using a 6-31G** basis set (Gaussian 98).40 The total nonbonded potential interaction energy of the simulated system consisted of long-range Columbic interactions between partial atomic charges calculated using the Ewald summation technique within an accuracy of 10-6 kcal/mol and van der Waals interactions computed using an atom based summation technique. The cutoff distance for the van der Waals interactions was kept at 9.5 Å. The bonding interactions were modeled by a composite force field to reflect the hybrid nature of the organic-inorganic system. The inorganic part is not dynamic and requires a hard potential, while the organic component requires a softer potential that can reflect its dynamic nature. The potential energy was computed using the universal force field (UFF)41 for the PbI4 layers, while the alkylammonium chains modeled using the polymer consistent force field (PCFF).42 PCFF parameters are optimized for predicting the structure and conformation of alkyl chains. Two nanosecond MD calculations were performed with the Material Studio package using the Forcite Plus module, which uses the velocity-Verlet integrator method for computing the positions and velocities of atoms.38 Simulations were initiated on a constant composition isothermal-isobaric (NPT) ensemble with a time step of 1 fs. Equilibrium values of the lattice parameters were judged to have been reached when these quantities fluctuate around their average values that remain constant over time. Typically, equilibrium values were reached within the first 100-120 ps. Structural data such as the lattice parameters are derived from the last 150 ps of the 2 ns trajectories. During the simulation, the temperature of the ensemble was maintained by the Andersen thermostat with a collision ratio of 1.00.43 To maintain the pressure an Andersen barostat with a cell time constant of 1 ps was used.43 The equivalent hydrostatic pressure was set to 0.1 MPa. Periodic boundary conditions were applied in three dimensions so that the simulation cell is effectively repeated infinitely in each direction. MD simulations were performed on the anchored bilayer for a range of temperatures from 65 to 665 K in steps of 30 K (see Supporting Information). Once the temperature
As described in Methodology the simulations were carried out using an NPT ensemble wherein the basal spacing as well as the in-plane lattice parameters were free to vary. Most of the bulk properties such as the volume, the lattice parameters and the total energy converged within the first 100 ps after which they fluctuated about a mean value. The system showed similar convergence of bulk properties at all temperatures (see Supporting Information). At convergence the value of the basal spacing was found to be close to the experimental determined spacing of (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) for all alkyl chain lengths. For example the room temperature interlayer spacing of (CnH2n+1NH3)2PbI4 is 33.6 Å while that obtained from simulations was 31.73 Å (see also Supporting Information). The in-plane structure of PbI4 layer shows a small deformation from the reported crystal structure, probably due to the small size of the supercell compared to the macroscopic dimensions of a true crystal. The deformation is however small and is ignored in the rest of the discussion (see Supporting Information). Thus an equilibrium state of the system is successfully attained from an arbitrary initial structure using the NPT ensemble. A snapshot of the ensemble at the end of the simulation shows that the anchored alkyl chains are arranged as a tilted bilayer with a herringbone arrangement (Figure 1). The chains tilt relative to the surface normal at an almost constant angle of 55°, thus optimizing interchain van der Waals interactions. The observed tilt angle is also in agreement with the value calculated from the ratio of the cross-sectional area of an all-trans alkyl chain, Ac, to the available surface area per alkyl chain, As, in the anchored bilayers of (CnH2n+1NH3)2PbI4. The tilt angle so determined, cos-1 (Ac/As), has a value of 60.8°. This value of the tilt angle would satisfy the criteria for formation of a tilted bilayer arrived at from theoretical considerations and would also exhibit a reversible phase transition.16 Phase Transitions. The (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) compounds exhibit a first-order endothermic differential scanning calorimetry (DSC) transition accompanied by a dilation of the lattice along the interlayer normal.19 The transition temperatures are reported to increase linearly with chain length. To investigate the thermal behavior of the simulated system we have monitored the temperature variation of three independent parameters (i) the heat capacity at constant pressure (Cp) of the system, (ii) the interlayer lattice spacing, and (iii) ∆B, the distance-fluctuation criterion for melting introduced by Berry et al.44 Heat Capacity. The heat capacity at constant pressure, Cp, is an important thermodynamic quantity that characterizes a phase transition, particularly melting.45 For the NPT ensemble Cp is related to the fluctuations in the enthalpy, 〈δH2〉
Cp )
〈δH2〉NPT 2
RT
)
〈H2〉 - 〈H〉2 RT2
The temperature variation in the computed values of the heat capacity of (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18) are shown in Figure 3a. The computed values of Cp show a gradual increase with temperature and exhibit an endothermic peak with the peak temperature (Tm) varying linearly with temperature (Figure 3b); a behavior similar to that reported from the experimental DSC profiles for this system. The experimental DSC curves of
Melting of an Anchored Bilayer
J. Phys. Chem. C, Vol. 114, No. 10, 2010 4539
Figure 3. (a) Variation of the heat-capacity, Cp, for different values of the alkyl chain length in (CnH2n+1NH3)2PbI4 as a function of temperature. (b) Variation of the transition temperature, Tm, with chain length (n).
Figure 5. Plot of the Berry parameter, ∆B, as a function of the reduced temperature, T/Tm, for different values of the alkyl chain length in (CnH2n+1NH3)2PbI4. The calculations are restricted to the atoms of the alkyl chains in (CnH2n+1NH3)2PbI4.
Figure 4. Variation of the interlayer lattice parameter c as a function of temperature for different values of the alkyl chain length in (CnH2n+1NH3)2PbI4. The transition temperature from the heat capacity data is indicated.
(CnH2n+1NH3)2PbI4 (n ) 12, 14, 16, 18), however, show an additional peak at low temperatures that has been attributed to the rotation of the headgroup.19 This endotherm is not observed in the simulations. It may, however, be pointed out that such a feature is not observed in the DSC of the closely related (CH3NH3)(CH3(CH2)nNH3)2Pb2I7 (n ) 11, 13, 15, 17) system.20 Although the simulations are able to reproduce the linear variation of Tm with chain length the temperatures are quite different from the experimentally reported values. For the n ) 12 compound the experimental Tm is 353 K while that returned by the simulations is 255 K. Similarly for the n ) 16 and 18 compounds the experimental transition temperatures are 364 and 373 K, respectively, while the calculated values are 365 and 455 K. This discrepancy in melting temperatures is due to inadequacies in the atom-atom potentials used and could, in principle, be rectified by fine-tuning the force-field parameters but since the objective of the present study was to understand the molecularity of the transition this was not attempted. Lattice Spacing. The variation of the computed interlayer lattice spacing, c, as function of temperature for the (CnH2n+1NH3)2PbI4 (n ) 12, 14, 16 and 18) compounds is shown
in Figure 4. The interlayer spacing shows a sharp increase at a temperature that is identical to that of the endothermic peak in the computed heat capacity vs temperature profile (Figure 3). The lattice parameter remains essentially constant both above and below this temperature. The simulations are, therefore, able to reproduce the experimentally observed jump in the X-ray diffraction determined basal spacing at the temperature of the endotherm in the DSC.19 Berry Parameter. To demonstrate that the observed transitions in Cp and the lattice spacing in both experiment and simulations (Figures 3 and 4) correspond to a melting of the alkyl chains of the anchored bilayer we have calculated the distance-fluctuation criterion for melting, ∆B, introduced by Berry.44 The Berry parameter ∆B is the relative root-meansquare (rms) fluctuation in the interatomic separation and is defined as
∆B )
2 N(N - 1)
∑ i
