Molecular Mechanics Study of Nickel(II) Octaethylporphyrin Adsorbed on Graphite(0001)

The effects of adsorption on the graphite(0001) surface on the nonplanar distortions of nickel(II)octaethylporphyrin were studied by molecular mechanics (MM) approach. Using the Consistent Force Field (CFF) program with previously developed parameters for metalloporphyrins and supplemented to treat intermolecular interactions geometry optimizations were carried out for 43 conformations of 28 distinct conformers of nickel(II)octaethylporphyrin. The stable energy-minimized conformers were stereochemically characterized, analyzed by the Normal-coordinate Structural Decomposition (NSD) method and compared with the available theoretical and experimental data for the isolated nickel(II)octaethylporphyrin structures.


Introduction
The chemistry of metalloporphyrins continues to draw attention not only because of their biological function, but also due to their current and potential application in many different industries. During the past decade the interaction of porphyrins with carbon materials has become a subject of numerous studies focusing on different physicochemical aspects. Important and extensively explored topics were the supramolecular assemblies of porphyrins on graphite [1][2][3][4], or fullerene [5,6], with an expected impact in the area of light-emitting diodes, organic displays, thin-film transistors, (photo)catalysts, design of new efficient photosynthetic systems, data storage media, and photovoltaic and electrochemical devices [7][8][9][10]. Although some theoretical work has been reported to analyze the adsorption mechanism [11][12][13], the conformations of porphyrin molecules adsorbed on a substrate have been the subject of only a few studies, for example, in connection with the calculation of STM manipulation signal of metalloporphyrins adsorbed on a metal surface [14,15].
Multiple lines of evidence [16,17] suggest that the flexibility of porphyrin core and puckering modes are key factors that determine and control the functional properties of metalloporphyrins. The unperturbed porphyrin macrocycle is conjugated and therefore expected to be planar. Non-planar geometries arise as a consequence of the electronic and steric effects of peripheral substituents and the nature of metal-ion macrocycle interaction [18]. In addition, as a result of adsorption process, molecule-surface interaction may induce conformational changes of the porphyrin core. Despite the size and flexibility of the porphyrin macrocycle, only a limited number of out-of-plane (oop) deformations are frequently seen. For substituted porphyrins a non-planar distortion can be classified on the basis of irreducible representations of the D 4h point group of a planar porphyrin [19][20][21]. The five most commonly observed distortions, represented by displacements along only the lowest-frequency normal coordinates, are: ruffling (ruf, B 1u ), saddling (sad, B 2u ), doming (dom, A 2u ), waving (wav, E g ), and propellering (pro, A 1u ) [22].
The roles of different types of out-of-plane distortions in biologically relevant metalloporphyrins are commonly studied on stereochemically restrained models of metalloporphyrins that have chosen conformations similar to those observed in proteins. Our previous studies [23,24] of octa-and tetra-halogeno tetraphenylporphyrins and their Ni(II) and Tb(III) complexes have shown that the type and degree of nonplanar deformation can be controlled by the peripheral substitution pattern, the steric bulkiness of substituents and the size of central metal of the macrocycle. In this work, we extend our studies by focusing on the changes of the porphyrin core conformation upon adsorption of metalloporphyrins on a graphite layer. Here we report a molecular mechanics (MM) study of the influence of a graphite(0001) surface on the nonplanar distortions of Ni(II)octaethylporphyrin, Ni(OEP), adsorbed on it, as well as the comparative MM study of isolated Ni(OEP) conformers using the previously developed force field [23,24] supplemented with new function and parameters which describe intermolecular interactions between porphyrin macrocycle and the graphite(0001) layer. Ni(OEP) has been chosen due to its rich stereochemistry and diversity of rotamers, and due to the fact that it has been thoroughly investigated in the recent past, both theoretically and experimentally, so that Ni(OEP) may be considered to be a "reference molecule" for any methodological development in the computer modelling of metalloporphyrins.

Stereochemistry
Ni(OEP) has interesting conformational properties arising from different orientations of the eight ethyl (Et) groups. Neglecting the hindered three-fold rotation about the terminal C-C bond (sp 3 -sp 3 ) of the Et group, and assuming that C β -C Et rotation is essentially two-fold, it is possible to generate a total of 2 8 = 256 rotamers of Ni(OEP), 28 of which are unique nonredundant conformers (see Figure 1). They are classified on the basis of the number of Et groups pointing on one side of the mean porphyrin plane (arbitrarily labeled as α) or on the opposite side of this plane (labeled as β) into five classes: α 8 or β 8 , α 7 β 1 or α 1 β 7 , α 6 β 2 or α 2 β 6 , α 5 β 3 or α 3 β 5 , and α 4 β 4 , comprising 2, 16, 56, 112, and 70 rotamers, respectively, or 1, 1, 6, 7, and 13 distinct conformers, respectively. Figure 1 shows also the numbering convention for 28 distinct conformers (in bold, above the structural diagrams), as well as their degeneracies (in parentheses, below the structural diagrams).

28
(2) Although there has been a lot of structural [25][26][27] and spectroscopic investigations and theoretical work on Ni(OEP), to our knowledge this is the first detailed MM study of energies and geometries of all 28 unique conformers.
All the stable structures, obtained by MM calculations were stereochemically characterized and analyzed by the normal-coordinate structural decomposition (NSD) method [22]. The Ni(OEP) conformers adsorbed on graphite surface were compared to the corresponding gas-phase conformers of Ni(OEP), which were in turn compared to the available X-ray structures of Ni(OEP).

Intramolecular Potential
Molecular mechanics calculations were performed with the 2007/PC version of the Consistent Force Field (CFF) conformational program [28]. Conformational energy was defined in the usual way as: where the potential functions and parameters for the first four terms (representing summations over all energy contributions for the isolated nickel(II)octaethylporphyrin) are described previously [23,24]. In this study we used the same force field parameters and functions, with the exception of the non-bonded potential, which was treated with Lennard-Jones 12-6 function instead of 9-6 function used previously. Parameters of the present 12-6 function were least-square fitted to reproduce the same r * and ε * values. This modification was introduced in order to get a more balanced ratio of numerical values for the energy terms (Eqn. 1), in other words, to ensure more reasonable scaling of variables, which is known to improve the rate of approach to the minimum of the total energy [29] in geometry optimizations. This choice of the non-bonded function produced only insignificant changes in the results of our previous studies [23,24], and did not affect any of our previous conclusions.

Modelling of Graphite Layer
One layer of the graphite(0001) surface, located in the xy plane was built up as a rigid neutral polyaromatic hydrocarbon (C 932 H 84 ) rectangular mesh (approx. 46×49 Å). During the minimization calculation the positions of the graphite atoms were kept fixed but the porphyrin molecule was allowed to move freely in all degrees of freedom (three translations and three rotations) in addition to the full relaxation of its internal degrees of freedom. Non-bonded cutoff was treated with a cubic spline switching function with the spline-on distance of 7 Å and the spline width of 4 Å.
The E (intermol.) NB contribution to the total energy consisted therefore of the sum of van der Waals and electrostatics interactions between the porphyrin macrocycle and the graphite C atoms. Van der Waals interactions were modelled using the Lennard-Jones 12-6 potential function with the same atom-specific parameters and combination rules as for the intramolecular nonbonded interactions [23,24].
Intermolecular electrostatic interactions were treated as monopole-quadrupole interactions between point charges located on the atomic positions of all metalloporphyrin atoms and uniaxial quadrupoles defined on each C atom of the graphite surface. Such a model was adopted due to the failure of the fixed atom-centered point charge description of the graphite surface in our modelling experiments, and the fact [30][31][32][33][34] that each C atom in graphite has an effective quadrupole moment. In the present force field the quadrupoles on graphite C atoms were constructed by placing negative charges (−q C ) along the normal to the graphite surface at ±a Å of each C atom, counterbalanced with the atom-centered positive charge +2q C (see Figure 2). The values of 0.5 Å for a, and 0.5 a.c.u. for q C were employed. This charge distribution resembles the one used by Vinter [35] in his XED (extended electron distribution) force field description of the C atoms in benzene.
A similar description of atom-centered multipoles could have been applied to the porphyrin core atoms of Ni(OEP) which exhibit aromaticity. However, the primary aim of this work was to compare the conformations of free Ni(OEP), for which we developed and optimized a force field based on atomcentered point charges previously [23,24], with the conformations of Ni(OEP) adsorbed on the graphite surface. As a consequence, we adopted a hybrid approach with graphite C atoms treated as quadrupoles and metalloporphyrin atoms treated as point charges. The monopole-quadrupole interaction energy, E MQ i j was calculated using the Cartesian form of the equation adapted from Hirschfelder [36]: where q i is the point charge on the i-th Ni(OEP) atom, r i j is the interatomic distance, e is the unit vector along r i j (i.e., r = |r| and e = r/r, see Figure 2), Θ j is the quadrupole moment tensor of the j-th C atom on the graphite surface, with components Θ αβ = i q i α i β i , for α, β = {x, y, z}.  (0001) surface.

Geometry Optimization
All the stable structures for isolated ("gas-phase") 28 conformers of Ni(OEP) were obtained by energy minimization starting from the planar structure, as well as from the four idealized non-planar forms (sad, dom, wav, and ruf ), which represent the normal deformations of the porphyrin core. They were generated from standard bond lengths and angles, and the corresponding z-coordinate displacements. Rotation of ethyl groups around the pyrrole-Et bond to achieve more favorable orientations was not observed. It was therefore possible to optimize geometries of all 28 conformers individually.
For the porphyrin macrocycle adsorbed on the graphite surface, all the stable conformers were found by energy minimization starting from various initial structures for each of the 43 conformers of Ni(OEP). The initial configurations were generated considering: (i) five conformations of the porphyrin core -as above, (ii) two different positions (one with metal atom directly above a given graphite carbon atom, and the other with metal located above the hole of the graphite hexagon), (iii) two different orientations of the porphyrin macrocycle relative to the graphite plane (one with the M-N bond in the porphyrin core eclipsed with respect to the C 2 axis passing through the C atoms of the graphite hexagon, and the other orientation staggered with respect to the same two lines), and (iv) various intermolecular distances (range 3-10 Å, step 0.5 Å). Considering that the potential energy surface for the porphyrin-graphite interaction is expected to be rather flat and shallow, this choice of initial configurations ascertained a reliable spanning of the conformational space for the porphyrin-graphite adduct. As in the case of isolated Ni(OEP) structures, the rotation of ethyl groups around the pyrrole-Et bond was not generally observed: only if the starting geometry has adjacent Et groups in a highly strained unfavorable orientation Ni(OEP) relaxes through Et groups rotation.
Geometry optimizations were carried out using the combination of steepest-descent, Davidon-Fletcher-Powell and Newton-Raphson methods [23,24,28]. Steepest descent and Davidon-Fletcher-Powell methods were mostly used, in that order, for initial exploratory searches and minimizations of conformations far from equilibrium [23,24,28]. The number of iterations varied widely in optimization experiments. In particularly difficult cases it was necessary to alternate between these two procedures more then once. To approach true minima Newton-Raphson iterations were always employed. Geometry optimizations were carried down to the energy rms gradient of < 10 −6 kJ/molÅ.

Normal-coordinate Structural Decomposition
For each of the equilibrium structures obtained by the energy minimization procedure we have performed normal-coordinate structural decomposition (NSD) analysis using the procedure of Jentzen, Song and Shelnutt [22], and the software available at http://jasheln.unm.edu. In this procedure distortions of the 24 atoms of a porphyrin core from the ideal D 4h symmetry are very adequately described as distortions along the lowest-frequency normal coordinates.

Nickel(II)octaethylporphyrin, Ni(OEP) (isolated)
For all theoretically possible conformers the energy minimization and geometry optimization procedure resulted in a unique stable structure, which did not depend on the choice of the initial nonplanar deformation of the porphyrin core. Structural parameters for selected resultant equilibrium conformations, together with the corresponding crystallographic data [25][26][27] are presented in the Appendix Table 1A. The resultant equilibrium conformations were also compared to those of Stoll et al. [37] who performed DFT calculations on Ni(OEP) and its isotopomers, and this comparison is also given in Table 1A. As can be seen, conformers labelled 19, 18 and 16 correspond clearly to the Triclinic A, Triclinic B, and Tetragonal forms of Ni(OEP), respectively, and the geometry of the resultant conformers is in good agreement with the one reported in the X-ray crystal structures [25][26][27]. Relative energies (graphically depicted in Figure 3(a)), energy contributions, Calculated Boltzmann population, and the results of the normal-coordinate structural decomposition (NSD) for all 28 conformers are summarized in Table 1. Table 1. Relative minimum energies ∆E and energy contributions (in kcal/mol) from bond stretching (E b ), angle bending (E θ ), torsional (E φ ), van der Waals (E vdW ), and Coulomb (E c ) interactions for the 28 isolated Ni(OEP) conformers, calculated Boltzmann population (P) at 298K, and out-of-plane (Doop) distortion of the porphyrin core (in Å). Data for the global minimum (conformer 26) are italicized. As can be seen from Table 1 and Figure 3(a) each of the 28 equilibrium conformers of Ni(OEP) can be assigned to one of the five groups of conformers (A, B, C, D, E) with nonoverlapping energies, and characterized by the number of pyrrole rings with ethyl groups oriented in the same direction. The groups are defined (referring to Figure 1) as: A, with four pyrrole rings with the same orientation of ethyl groups (conformers 1, 8,16,19) and the energy in the range from −17.11 to −17.40 kcal/mol; B, with three such pyrrole rings (conformers 2,9,14,15) and with energy from −17.57 to −17.86 kcal/mol; C, with two such pyrrole rings (conformers 3,4,5,6,7,17,18,20,21,22,23,24) and with energy from −17.92 to −18.29 kcal/mol; D, with one pyrrole ring with the same orientation of ethyl groups (conformers 10, 11, 12 and 13) and energy from −18.35 to −18.63 kcal/mol. The lowest energy conformers, group E, (energies in the range −18.71 to −19.11) have opposite orientation of the ethyl groups (conformers 25, 26, 27 and 28) on all pyrrole rings. The mutual orientation of ethyl groups on neighboring pyrrole rings does not have any significant influence on the energy value, but does influence the total core puckering.  Although the total energies between all conformers are slightly different, different orientation of ethyl groups cause a distinction in the degree of non-planarity and NSD pattern (Table 2A). Thus, conformer labelled 26 has the greatest puckering amplitude (pure ruffled conformation), and this conformer is the global minimum for Ni(OEP) species, while the conformer 28 is the most planar one. This is in agreement with a DFT calculation [37], in which the ruffled conformation was shown to be energetically favored for more than 0.2 kcal/mol. When the substituents on the porphyrin have more than one possible combination of orientations, then the conformation that occurs in the crystal depends on the relative energies of the conformers. If energy differences among different stable conformers are large, the conformation observed in the crystal is likely to be the most stable one. If the energy differences are small, several conformations may be observed in the crystalline state. The relative energies of the stable conformers, obtained by the present MM calculations, indicate that all considered conformers may appear in crystalline state as well as in the solution.

Nickel(II)octaethylporphyrin Adsorbed on Graphite
Results are presented for the 43 conformers of Ni(OEP) adsorbed on the graphite surface. Selected resultant equilibrium conformations are shown in Figure 4. Calculated minimum energies are given in Table  2, and graphically depicted in Figure 3(b). Selected NSD results compared to the NSD results of corresponding isolated calculated structures are shown in Figure 5. Complete NSD results for all Ni(OEP) conformers adsorbed on the graphite surface, as well as the results for isolated Ni(OEP) conformers are summarized in Appendix Table 2A. First fifteen conformers (Figure 1), differ in number of ethyl groups lying above(α) and below (β) porphyrin mean plane, were considered in two orientations: when α or β ethyl groups are pointed toward to the graphite surface. Since the others conformers have the equal number of ethyl groups on the both side of porphyrin plane we consider only one orientation of ethyl groups with respect to the graphite surface. MM calculations resulted in a unique stable conformation for all conformers, regardless of the initial nonplanar deformation of the porphyrin core, the initial relative orientation of porphyrin macrocycle and graphite layer, and their initial distance. All stable conformations obtained for 43 conformers differ in energy (Table 2 and Figure 3(b)), in core puckering ( Figure 5, Table 2A) and in the position of porphyrin relative to graphite. As can be seen from Table 2 significant difference in the total energy value is due to the contribution of intermolecular quadrupole-monopole interaction. The analysis of the quadrupolemonopole contributions to the total strain shows that the energy linearly decreases, with the increasing number of ethyl groups pointed toward to the graphite surface. Thus, the lowest energy structure is α 8 , while β 8 has the greatest energy value. Figure 3(b) shows 9 distinct energy groups which differ in the number of ethyl groups pointing towards the graphite surface. The overlap of some groups is due to the fact that the total energy (a sum of intra-and intermolecular interactions) is presented.
Possible interactions between porphyrin molecule and the graphite layer are π-π (π-stacking), σ-π (Et-graphite interactions that can become repulsive at small distances), and M-π. The conformer with all Et groups opposite to the graphite layer (α 8 ) has dominantly π-π interactions with graphite; the one with all Et groups pointing towards the graphite surface (β 8 ) has σ-π interactions; the others possess combinations of both.
In comparison to the isolated structures, porphyrin cores are more puckered for all conformers, except for α 8 adsorbed on graphite ( Figure 5, and Table 2A), with the presence of dom deformation. This enhanced puckering is a result of the balance between π-stacking interactions that tend to flatten the porphyrin core, and repulsive forces involving interactions between Et group or the central metal atom with the graphite C atoms. The conformer α 8 adsorbed on graphite layer is less puckered, however, presumably due to the fact that all ethyl groups point away from the graphite surface enabling the porphyrin core to approach the surface more closely and to enhance the flattening π-stacking interactions. Another consequence of adsorption of conformer α 8 on the graphite surface is the increased doming of the porphyrin core due to Ni-C(graphite) repulsions, which displaces the Ni atom away from the graphite surface. At distances greater than 9.5 Å there are no appreciable intermolecular interactions and the NSD patterns are the same as for the isolated porphyrin molecule. In all simulations a rotation of ethyl groups towards more favourable orientations was not observed.

Movement of Ni(OEP) in the vicinity of graphite(0001) surface
The primary driving force for adsorption of a particular Ni(OEP) molecule onto the graphite(0001) surface in the present MM modelling is the long-range dispersion interaction between the two moieties. Geometry optimization always leads to the parallel (π-π stacked, or "face-to-face") orientation of Ni(OEP) on graphite. However, contrary to what might be assumed, Ni(OEP) molecules do not approach the surface in an unvaryingly parallel manner even when an optimization starts from a perfectly parallel orientation. The approach of Ni(OEP) is rather characterized by swinging or rocking movements of Ni(OEP) as various atoms on the periphery of the porphyrin core approach the graphite layer, eventually ending in a parallel orientation (Figure 4).
In addition, the approach of Ni(OEP) is accompanied with lateral movement, which is small in accordance with the dissimilarity between the size of the unit pattern of the graphite surface and that of the Ni(OEP) molecule.
In view of the fact that benzene molecules can adopt other than parallel mutual orientations [38] we also performed geometry optimizations starting with a T-shaped constellation of the two moieties. Very slow initial convergence (due to a reduced number of long-range interactions between Ni(OEP) and graphite atoms) eventually led to the parallel orientation through a rolling of the whole Ni(OEP) molecule.
Neither the lateral movement nor the rolling and turning of Ni(OEP) was possible in a force field based on simple monopole-monopole electrostatic interactions between the moieties, which corroborated the need to model graphite C atoms as axial quadrupoles.

Concluding remarks
The adsorption of Ni(OEP) species on graphite(0001) layer was analyzed as model case between porphyrin molecules and chemically inert surface. The intermolecular interactions were modelled using the Lennard-Jones 12-6 potential functions and monopole-quadrupole electrostatic interactions.
We have shown that adsorption on a surface is an additional factor that should be taken into account in conformational analysis of metalloporphyrins. MM calculations and NSD analysis revealed that isolated Ni(OEP) structures, and Ni(OEP) structures adsorbed on the graphite layer, differ in core puckering.
It is well-known [16] that the type and magnitude of normal deformations has profound consequences on spectral, electrochemical and other properties of porphyrins. Thus, changes in physical and chemical properties, as well as metalloporphyrin functionality, when it is adsorbed on a surface is a consequence not only of adsorption (and the presence of the surface), but also of specific conformational changes.
Scudiero et al. [39] in their scanning tunnelling microscopy (STM) investigation of Ni(OEP) on a highly ordered pyrolytic graphite found that Ni(OEP) self-assembles on the graphite surface in the form of a flat 2D lattice. In agreement with this experiment we determined that parallel mutual orientations are always favored irrespective of initial orientation of Ni(OEP). Since STM technique cannot directly show the orientations of ethyl groups, the present MM approach might be a useful complement in structure determination and in the elucidation of self-organization of porphyrins on solid substrates.

Acknowledgments
This work was financially supported by the Serbian Ministry for Science and Environmental Protection through the Grant No. 142017G. We thank Professor Carlo Adamo (ENSCP, Paris) for the initial stimulus as well as for the constant interest and advice.