Molecular simulations for the conformational assessment of a porphyrin–fullerene dyad in different environments

Conformational space of a porphyrin–fullerene dyad with the donor and acceptor connected by a relatively flexible linker is studied by molecular dynamics simulations in both non-polar and polar solvents, as well as in vacuum.

The most probable conformations obtained from the vacuum MD simulations were optimized with semi-empirical (SE) and density functional theory (DFT) methods and the extent of the structural changes is assessed.

The computational results indicate the co-existence of different conformers in both polar and nonpolar solvents showing agreement with experimental results.

The most probable vacuum conformations at 300 K are similar to the ones at 0 K, while the structures most often observed in the solvents show less compact conformations.

Optimization with SE and DFT calculations leads to structures, which represent relatively well the folded conformations in solvent, which validates the electronic structure calculations relevant to describing photoinduced electron-transfer in H2P–O34–C60.


In recent years, photoinduced electron-transfer (ET) processes taking place in covalently linked porphyrin–fullerene dyads have attracted special attention.1–4

In such systems photoexcitation of porphyrin (donor) leads to a chain of ET processes ending at the fullerene (acceptor).

The spatial arrangement of donor–acceptor (D–A) molecules is a subject of primary interest, as the rate of ET depends considerably on the intermolecular distance and relative orientation of these molecules.

In the present study we are searching for the possible conformations of a porphyrin–fullerene dyad (see Fig. 1) in which the D–A molecules are connected by a single relatively flexible linker.

This allows for a substantial degree of conformational flexibility in the molecular topology and two limiting scenarios evolve: (a) a folded conformer in which the fullerene molecule is oriented on top of the porphyrin molecule, and is therefore, characterized by a small D–A distance or (b) an extended conformer in which the fullerene molecule is oriented away from porphyrin, and is therefore, characterized by a relatively large D–A distance.

Spectroscopic measurements indicate the co-existence of different conformers in both polar and nonpolar solvents and a two conformer model was proposed to account for the experimental observations.3

In addition to the spectroscopic observations, there is evidence from NMR measurements performed in chloroform on the Zn analog of the dyad in this study, which indicates that the folded conformer is the dominating one or that several conformations are rapidly interchanging.5

In order to study the conformational behavior of the porphyrin–fullerene dyad in different environments, molecular dynamics (MD) simulations were performed both in non-polar (C6H6, benzene, dielectric constant, εr = 2.28) and polar solvents (C2H6SO, dimethyl sulfoxide, DMSO, εr = 47.2).

In addition, MD simulations were performed in vacuum with the purpose of assessing the differences between solvent and vacuum structures, and therefore, the usability of the vacuum structure for the more detailed electronic structure calculations relevant to the photoinduced ET in H2P–O34–C60.

The most probable conformations obtained from the vacuum MD simulations were optimized with semi-empirical (SE) and density functional theory (DFT) methods and the extent of the structural changes was assessed.

Computational methods

Molecular dynamics simulations

The GROMACS (version 3.1.4) molecular dynamics simulation package with the GROMACS force field (based on the GROMOS-87 force field parameters6,7 with some modifications8,9) was used for the studies.

Periodic boundary conditions were applied to minimize the artificial effects caused by the finite size simulation box (approx.

70–110 nm3).

Atomic detail was used except for hydrogen atoms bound to the carbon atoms of the porphyrin–fullerene dyad.

They were treated as united atoms.

In addition to the 146 atoms of the dyad, the simulation box contained 530–923 solvent molecules resulting to the total of about 2300–3800 (in DMSO) or 6500 (in benzene) atoms depending the solvent and starting structure of the dyad.

The simulations were performed with a 2 fs time step with all bond lengths constrained with the LINCS10 algorithm.

The solvents were included explicitly in the calculations to make the environment as realistic as possible for the molecule of interest.

The particle-mesh Ewald11–12 method was applied to the electrostatic interactions for the simulations to avoid the artefacts created when cutting off the long range electrostatic interactions.

The charge group-based pair list was updated every ten steps.

The temperature was controlled using weak coupling to a heat bath of 300 K with a time constant of 0.1 ps.13

Similar procedure with a time constant of 1.0 ps was applied to keep the pressure at 1 bar.

Initial velocities of the atoms were randomly generated from the Maxwell distribution at 300 K.

First, after the solvation of the porphyrin–fullerene dyad in the case of DMSO and benzene simulations, the system was energy-minimized (EM) with the steepest descent method until the maximum of the force gradient components was smaller than 100 kJ mol−1 nm−1 or the maximum of 1000 steps was reached.

Then, the solvent was equilibrated at least for 10 ps while restraining the dyad using harmonic potentials with the force constant of 1000 kJ mol−1 nm−2 to remove solvent based bad contacts.

After the equilibration of the solvent, the production runs were performed for 50 ns.


The semi-empirical Austin model (AM1)14 and the parametric method number 3 (PM3)15,16 were employed for geometry optimization.

Both parametrizations are known to provide reliable ground-state geometries for porphyrins and porphyrin–fullerene dyads.17–20

The AM1 and PM3 optimized geometries and the corresponding energies were calculated with Gaussian.21

The AM1 and PM3 single point energies for the two MD conformers were computed with AMPAC software22.


The Perdew–Burke–Ernzerhof (PBE)23–26 generalized-gradient approximation (GGA) exchange-correlation functional was employed for geometry optimization.

The Karlsruhe split-valence basis set augmented with polarization functions (SVP)27,28 was used.

The SVP basis set consists of two basis functions for H (4s)/[2s], six basis functions for C (7s4p1d)/[3s2p1d], and six basis functions for N (7s4p1d)/[3s2p1d].

The terms in parentheses represent the numbers of primitive functions of each type and the terms in square brackets represent the numbers of contracted basis functions of each type.

The resolution of identity approach (RI) was used to reduce the computational effort, and accordingly, the time needed for calculations.29

All calculations were performed with the Turbomole software package30.

Results and discussions

Molecular structures: MD study

In vacuum

Fig. 2(a) shows the distributions of different atomic and atom group distances within the porphyrin–fullerene dyad obtained from 50 ns MD runs in vacuum at 300 K applied to two different starting conformations: a folded [Fig 2(b)] and an extended one [Fig. 2(c)].

Similar distributions [Fig. 2(a)] were generated irrespective of the starting conformation.

Between the centres of the fullerene and porphyrin only one optimal distance of about 0.65 nm can be deduced (one peak in the distribution), while the linker tends to show two favourable conformations (two peaks or a peak with a shoulder in the corresponding atomic distance distributions: N20–C84, N64–C84, N20–C82, N64–C82).

In Fig. 3(a) and (b), the two peak conformations of the distributions of Fig. 2(a) are displayed: 3(a) corresponding to the higher and 3(b) to the lower peaks of the atomic (N–C) distance distributions.

Thus, the double peaks of the N–C distributions indicate that the linker can bend either from left (close to N20) or from right (close to N64) over the porphyrin.

The porphyrin ring has a ruffled conformation in both cases.

In a ruffled conformation the meso carbon atoms are alternately above and below the porphyrin mean plane (least-squares plane calculated for the carbon and nitrogen atoms of the ring) while the core nitrogens are in the plane.31

The ring inversion of the cyclohexane ring in the linker, which is mainly responsible for the conformations presented in Fig. 3(a) and 3(b), has been previously observed experimentally in 1H-NMR measurements of a ZnP–O34–C60 compound,5 that is a Zn analog of the H2P–O34–C60 molecule.

Fig. 3(d) depicts the frames obtained every 500 ps from the trajectory of a 50 ns MD run.

The conformations obtained in vacuum at 0 K look similar to the peak conformations of the 300 K distributions shown in Fig. 3(a) and (b).

If any of the vacuum conformations appearing at 300 K is taken as a starting structure for a simulated annealing (SA) from 300 to 0 K, the SA run leads to a structure similar to one of the two structures shown in Fig. 3(a) and (b).

In polar solvent

The distributions of atomic and atom group distances within the porphyrin–fullerene dyad obtained from the 50-ns MD runs in DMSO are shown in Fig. 4.

Both folded and extended starting structures, similar to the ones shown in Fig. 2(b) and (c), were considered.

Fig. 4. shows distinct starting structure dependent differences in the conformational distributions indicating that a single 50 ns run does not exhaustively sample the conformational space.

Repeating runs with different initial conformations enhances the conformational search and make the sampling of the large conformational space more feasible for the available computing capacity.

It appears that there are at least three energy minima in the conformational energy hypersurface in DMSO at 300 K (Fig. 5), two close to the folded starting structure [Fig. 4(a)] and a broader one with a slightly more open conformation [Fig. 4(b)].

The structures in Fig. 5 were obtained by applying a procedure which first selects a set of conformations with the porphyrin–fullerene distance at the peak of the distribution (the solid curves in Fig. 4) and then prunes the set until only those structures which also correspond to the peaks of the distributions of all the curves in Fig. 4 were left, i.e. in the selected structures the linker conformations were also the ones most often detected during the simulations.

The distance (i.e. the position of the peak of the distribution) between the centre of the fullerene and the centre of the four nitrogens of the porphyrin ring is 0.87 nm for the higher peak and 0.68 nm for the lower peak of Fig. 4(a).

In the third structure [Fig. 4(b)] the corresponding fullerene-porphyrin distance is 1.31 nm (and 1.19 nm for the shoulder of the peak).

The average distances are also shown in the figure.

The main difference between the conformations of Fig. 4(a) and (b) (or between the light grey structure and the other ones in Fig. 5) is the 180° rotation of the torsional angle around the OC-phenyl bond (see Fig. 1).

The rotation explains the notably different peak positions in Fig. 4(a) and (b), while the broad or double peaks within one curve are a result of the sidewise movement of the fullerene and/or the linker on the porphyrin.

Torsional rotations were not detected around the other single bonds of the linker (N-phenyl and phenyl porphyrin) during any of the simulations.

In non-polar solvent

The MD results of the 50 ns runs in benzene are shown in Figs. 6 and 7.

As before, simulations starting from both folded and extended conformations were performed.

Fig. 6 depicts the distributions of the atomic and atom group distances within the porphyrin–fullerene dyad during the simulations.

The centre to centre distances at the top of the peaks (solid curves) were 0.89 nm and 1.37 nm for Fig. 6(a) and (b), respectively.

The structures in Fig. 7 were obtained by applying the similar procedure to the distributions of Fig. 6, as the one used to find the structures of Fig. 5 from the distributions of Fig. 4.

Fig. 7 illustrates that the top of the peak of the centre to centre distribution of Fig. 6(b) includes significantly different conformations.

Very similar conclusions to the ones presented above for Fig. 4 can be drawn from Fig. 6: The main differences between the peak position of Fig. 6(a) and (b) are related to the torsional rotation around the OC–phenyl bond of the linker, while the double (or broad) peaks are a result of the sideways motions of the fullerene and/or linker on the porphyrin.

Solvent effects

In addition to the MD simulations at 300 K, some simulated annealing runs were also performed from 300 to 0 K in the solvents.

It appeared that although the starting structures of the SA runs differed only slightly at 300 K (e.g. highly similar frames very close to the top of the peaks of the distributions of Figs. 4 and 6), they hardly ever ended up to the same structure at 0 K after SA.

In fact, the SA simulations in solvent usually led the dyad to highly variable conformations, even though the starting structures were very close to each other.

Only the porphyrin ring always maintained the ruffled conformation.

Thus, it can be concluded that the energy hypersurface in the solvents appears to be rough at the bottom, i.e. at 0 K, while only a few (broader) local energy minima were found at 300 K. However, quantitative energetic evaluations could not be performed for the different energy minima of the dyad in solvent due to the massive background noise generated by the solvent–solvent interactions of the system.

The shorter atomic distances obtained from the vacuum MD simulations are related to missing interactions of the dyad with the solvent molecules.

The fact that vacuum simulations tend to produce compact globular conformations have frequently been discovered in different systems previously.32,33

The significantly larger centre to centre porphyrin–fullerene distances obtained for the simulations in the solvents starting from the extended conformation suggest that the solvent molecules have relatively strong interactions with the dyad.

The formation of fully folded conformations similar to the ones observed in the vacuum simulations (Fig. 3) has steric hindrance caused by the solvent molecules surrounding the dyad, and thus, a partial desolvation is needed for the dyad to attain a compact structure like the ones in Fig. 3.

Furthermore, complete sampling of the conformational space is not guaranteed to be reached with the present 50 ns simulations, but would be reached with a long enough run.

Especially in the case of benzene, a relatively viscous liquid, sampling in a single MD run is likely to be more or less confined near to a few local minima, with the exact region of the conformational space sampled depending on the starting structures and velocities.

Molecular structures from SE and DFT compared with the MD results

The 0 K vacuum MD conformations (obtained by performing, first, a SA run from a 300 K structure to 0 K, and then, an EM) similar to the ones shown in Fig. 3(a) and 3(b) were optimized with SE methods using AM1 and PM3.

In addition, the conformation in Fig. 3(a) was also optimized using DFT.

The two distinct conformations of the linker which are mainly differentiating the starting structures (see Fig. 3(a) and 3(b)) are retained also after optimization at SE level and will most probably be retained also after DFT optimization or other higher level of calculation.

In all cases the porphyrin ring remained in a ruffled conformation, the degree of ruffling being least in the case of the structures optimized at SE level and highest for the representative MD structures.

The mean deviation of the porphyrin ring atoms from the least-squares plane calculated for the carbon and nitrogen atoms of the ring is 0.02 nm for the representative MD structures and just about half that large in the case of the DFT optimized structure.

The SE optimizations yield structures in which the porphyrin ring is more planar with the mean deviation only about 0.005 nm.

The distances between the centre of fullerene and the centre of the porphyrin ring, as well as the distances between the atoms N64 and C82, N20 and C82, N64 and C84 and N20 and C84, are given for the two conformers in Table 1.

In addition, a graphical representation of these distances together with the results obtained form the different MD simulations is shown in Fig. 8.

Inspection of Fig. 8 and Table 1 reveals that independent of the computational method the differences between the centre to centre distance of fullerene and porphyrin in the two conformers are smaller than 0.01 nm.

However, the centre to centre distance differs by more than 0.15 nm between the structures optimized at different levels of approximation.

Fig. 8 shows that the centre to centre as well as the atom pair-distances of the MD simulations in the solvents starting from the folded conformation are relatively close to the corresponding SE results.

Only the vacuum MD results produced somewhat smaller centre to centre distances which, however, are not far from the DFT results.

Based on the data given in Fig. 8 and Table 1 we can conclude that the distance between the centre of fullerene and the centre of porphyrin is decreasing depending on the computational method used in the following order AM1 > PM3 > DFT > MDvacuum.

After SE optimization the fullerene molecule is located further away from porphyrin and the porphyrin ring is becoming more planar (see Fig. 8 and Table 1).

The centre to centre distance calculated with AM1 is larger by about 0.06 nm than the distance calculated with PM3.

The DFT optimization of the conformer in Fig. 3(a) yields a centre to centre distance which lies in between the one obtained from MDvacuum and PM3 calculations.

The same dependence on the computational method and decreasing trend is observed also for the distances between the atoms N64 and C82, N20 and C82, N64 and C84, and N20 and C84.

Based on the geometry optimization of the two conformers with Gaussian AM1, the difference in energy (heat of formation) between the two conformers is 0.42 kJ mol−1, which is insignificant.

The difference between the single point total energies computed with AMPAC (AM1 and PM3) for the MD optimized (0 K) conformers is ca.

55 kJ mol−1.

However, when the slightly different vacuum structures of 300 K were run to 0 K using, first, the SA run, and then, an EM, the porphyrin, linker and fullerene resulted to the similar structure, but the ligands of the porphyrin always showed some slightly different orientations.

The force field based energies of the vacuum MD structures at 0 K differed some 300–400 kJ mol−1 due to the small differences only in the ligands.

As this applies to both conformer 1 and 2, the possible energy differences between the two MD conformations (appearing smaller than the energy differences generated by the ligands only) cannot be distinguished, and thus, also the value of 55 kJ mol−1 obtained with AMPAC for the MD optimized conformers can be considered insignificant.


The MD studies of the porphyrin–fullerene dyad in vacuum and in two different solvents revealed that the most probable vacuum conformations at 300 K are similar to the ones at 0 K, while the structures most often observed in the solvents showed less compact conformations.

Both folded and more open, or extended, conformations were discerned from the solvent simulations the predominant structure depending on the initial conditions of the runs, as well as on the solvent.

When an extended structure was solvated as a starting structure for the simulations, the partial desolvation needed for the production of fully folded conformers appeared difficult due to the hindering solute–solvent interactions, and therefore, resulted in the dominance of more or less extended structures in those simulations.

Accordingly, the computational results are in agreement with the spectroscopic measurements3 indicating the co-existence of different conformers in both polar and nonpolar solvents and also validate roughly the suggested “two conformer” model of folded and extended structures.

The distance between the centre of the fullerene and the centre of the porphyrin ring in vacuum at 0 K is decreasing in the order of AM1 > PM3 > DFT > MDvacuum depending on the computational method used.

However, the SE (and DFT) optimized vacuum MD structures appear to represent relatively well the folded conformations in solvent, and thus, the electronic structure calculations based on SE (or DFT) can be expected to give reasonably realistic representation of other properties and phenomena including the photoinduced ET in H2P–O34–C60.