Quantum chemical replica-exchange umbrella sampling molecular dynamics simulations reveal the formation mechanism of iron phthalocyanine from iron and phthalonitrile

One-pot synthetic reaction mechanisms are desirable for large-scale industrial synthesis of chemical materials. Phthalocyanine (Pc) and phthalocyanine metal-complexes (MPcs), shown in Figs. 1(a) and 1(b), respectively, are examples for chemical materials utilized as pigments for decades due to their vivid green or blue color. For example, the copper phthalocyanine complex has a brilliant blue color, making CuPc an ideal photostable dye in pigments.¹ The Pc molecules and complexes possess a cyclic structure, in which all atoms are located in the same plane and are structurally similar to their cousins, the porphyrin molecules and complexes. Remarkably, MPc can be synthesized in a one-pot synthetic reaction, the so-called phthalonitrile (PN) method, albeit only with moderate yield.^2–4 It is a well-known synthesis in which Pc or MPc is obtained from a mixture of PN and metal powder at high temperatures above 280 °C.

Recently, MPcs have been used not only for dyes or pigments but also as organic semiconductors in covalent organic frameworks,^5,6 optical data storing media, etc. Somewhat surprisingly, Pc and MPc were discovered accidently in the early 20th century,^7–9 and even more surprisingly, the reaction mechanism underlying the nearly century-old PN method is not fully understood to date.¹⁰ 

Examples of hypothetical formation mechanisms for MPc from simple reactants in two different synthetic methods are shown in Fig. 2. In the Wyler method, which has been mainly used industrially due to its higher reaction yield, MPc is synthesized from phthalic anhydride with metal powder or metal salt and urea. The second is the phthalonitrile (PN) method, which is used on laboratory scale. In these methods, PN and metal powder or metal salt are heated up to 280 °C.^2–4 

To this day, the formation mechanism of Pc and MPc complexes has not clearly been established. Moreover, intermediates, transition states, and free energy barriers have not been identified, neither experimentally nor theoretically.

In order to theoretically study the mechanism for forming these molecules by the PN method, we performed quantum chemical molecular dynamics (QM/MD) simulations for the formation of iron phthalocyanine (FePc) from four PNs and one iron atom in the gas phase. Straightforward QM/MD simulations including periodic boundary conditions (PBC) on computationally manageable time scales (nanoseconds) do not yield the desired product, no matter how high a temperature is chosen. On the MD time scale, FePc formation is a rare event because high barriers on the free energy pathway, associated with the anionic cyclic polymerization of the four PN units, must exist. Unfortunately, microsecond time scales are not affordable even for approximate density functional theory (DFT)-based first principles MD simulations. To overcome this difficulty, we employed one of the generalized-ensemble algorithms in this work (for reviews, see, e.g., Refs. 11–15), which greatly enhance conformational sampling. These methods have traditionally been used in combination with classical Monte Carlo (MC) and molecular dynamics (MD) simulations. However, for the study of chemical reactions, some groups reported promising results from QM-MC and QM-MD simulations.^16–18 In the present study, we employed the replica-exchange umbrella sampling (REUS) method,¹⁹ which combines the replica-exchange method (REM)^20,21 and the umbrella sampling (US) method,²² in combination with the density-functional tight-binding (DFTB)²³ quantum chemical potential. This enhanced reactive MD method has been implemented²⁴ recently in the freely available DFTB+ program package.²⁵ The use of multiple umbrella potentials via the REUS method along the reaction coordinate is crucial for the study of the FePc self-assembly mechanism: Since it cannot be expected that the four PN molecules arrange themselves in a concerted single-step reaction around the iron cation, a stepwise bond formation dynamics with intermediate steps demands that multiple regions on the potential energy surface can be visited effectively. For this reason, the REUS with its enhanced sampling capability is required to reveal the intermediate steps of the self-assembly reaction.

This article is organized as follows. In Sec. II, the relevant methods are explained. In Sec. III, computational details of the simulations are given. In Sec. IV, the results of the simulations are presented. Section V presents the conclusions.

In the US method,²² we use an umbrella potential V(q), which is defined by

where ξ(q) is a reaction coordinate. Two adjustable parameters, a force constant k and a reaction coordinate midpoint value d, are the adjustable parameters in the umbrella potential. The umbrella potential V(q) is added to the original potential energy E₀, and the total potential energy is defined by

The umbrella sampling simulation allows the system to explore the region around ξ(q) = d, which may not be reached during regular MC or MD simulations due to the existence of high energy barriers.

In a conventional US simulation, M independent US simulations at a temperature T₀ are performed with M different umbrella potentials V_m (m = 1, 2, …, M). V_m (m = 1, 2, …, M) have midpoint distances d_m (m = 1, 2, …, M) and force constants k_m (m = 1, 2, …, M) which can depend on m and are defined by

Therefore, the total potential energy is given by

After M US simulations with M different E_m, the free energy and physical quantities of the original unbiased system, which is based on potential energy E₀ at temperature T₀, are calculated by using the weighted histogram analysis method (WHAM).²⁶ By solving the following WHAM equations self-consistently, the probability distribution P_T0(ξ) and the “dimensionless” Helmholtz free energy f_m of the original unbiased system can be obtained:^27,28

where N_m(ξ) is the histogram of distributions of reaction coordinate ξ and n_m is the total number of samples for the simulation with E_m. The potential of mean force (PMF), F_T0(ξ), which is the Helmholtz free energy as a function of the reaction coordinate of the original unbiased system, is given by

The expectation value of a physical quantity A of the original unbiased system is given by

The Replica-Exchange Umbrella Sampling (REUS)¹⁹ is one of the generalized-ensemble algorithms that combines REM and US. In REUS, the system consists of M non-interacting replicas i at M different umbrella potentials U_m(q) (i, m = 1, 2, …, M). The total potential energy E_m(q) by the mth umbrella potential is defined by Eq. (4), where

Here, the parameters $λmll=1,2,…, L; m=1,2,…, M$⁠, usually chosen empirically from short test simulations, are the weighting factors for the corresponding umbrella potentials. Their values are selected to ensure sufficient population overlap between those of neighboring umbrella potentials.

The sets of labels for replicas, i, and umbrella potential, m, are related by permutation,

where f(m) is a permutation function of umbrella potential label m and f⁻¹(m) is its inverse.

The state X in this generalized ensemble can be written as

where $xmi$ is specified by the coordinates $q[i]$ and momenta $p[i]$ of N atoms in replica i at temperature T_m,

Independent MC or MD simulation with M replicas is carried out simultaneously in the canonical ensemble at temperature T₀, and every few time steps, replicas i and j with neighboring umbrella potentials U_m(q) and U_m+1(q), respectively, are exchanged:

The replica exchange is performed by the following probability ω of the Metropolis criterion:

where

In a REUS simulation, by replica exchange, a random walk in reaction coordinate ξ is performed without getting trapped in local-minimum states. After a REUS simulation is performed, the PMF can be obtained as in US using the WHAM [see Eqs. (5)–(7)].

We employed the density-functional tight-binding (DFTB)-based QM/MD simulation with REUS, which is referred to as DFTB-REUS-MD. All calculations were performed using our modified version of the DFTB+ program.^25,24 The initial system was composed of four PN molecules and a single iron atom in gas phase, with initial coordinates arranged to resemble somewhat the arrangement of the molecular units in an FePc complex. The MD simulation time in these simulations was 150 ps per replica. To study the dependency of the trajectories on the initial structures, we also carried out a DFTB-REUS-MD simulation of 50 ps per replica by using different initial coordinates where PN molecules and the iron atom were placed more randomly as compared with the former simulation.

Regarding the DFTB, we employed the second-order self-consistent charge (SCC-DFTB, nowadays also referred to as DFTB2) flavor,²³ which can be considered as an approximate Density Functional Theory (DFT) method^29,30 within the Slater-Koster approximation.³¹ The SCC-DFTB method has reasonable accuracy and very low computational cost as compared to traditional DFT-based QM methods. The mio-1-1²³ and trans3d-0-1³² parameter sets were employed throughout this work. A finite electronic temperature of 1000 K was employed in order to obtain convergence of the DFTB electronic density, and the trajectories were therefore run on the DFTB Mermin free energy³³ surfaces.

The simulation of chemical reactions in a DFTB-REUS-MD trajectory is naturally sensitive to the initial placements of the reactant molecules. We therefore decided to perform two sets of simulations. In the first “A set,” we prepared four different sets of initial coordinates, which had four PNs near the iron atom and similar to FePc as shown in Fig. 3. In the second, “B set,” we placed four PNs around one iron atom randomly as shown in Fig. 4, where all replicas started from the same initial coordinates. After the simulation, we performed DFTB geometry optimization for several structures to quench the high-temperature structures.

During the formation of FePc, there were two types of covalent bonds between a carbon atom and a nitrogen atom. One is a C—N bond in “internal” PNs as shown in Fig. 5(a). The other is a C—N bond in “external” PNs, which means that two PNs are bound by a C—N bond as shown in Fig. 5(b). Therefore, we defined a reaction coordinate as a linear combination of eight distances between a carbon atom and a nitrogen atom in internal and external PNs as

Subscripts in Eq. (19) correspond to labels of atoms as shown in Fig. 6.

We used M = 16 replicas for the DFTB-REUS-MD simulations, which have different k and d values in each umbrella potential. Each replica had the same initial coordinates. In the “A set” DFTB-REUS-MD simulation, we set k values to be 1.0, 50.0, 100.0, 100.0, 100.0, 75.0, 50.0, 50.0, 50.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, and 20.0 kcal/mol Ų and d_m values were chosen to be 11.0, 11.5, 12.0, 12.25, 12.5, 12.75, 13.0, 13.25, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5, and 17.0 Å. In the “B set” DFTB-REUS-MD simulation, we set k = 1.5 kcal/mol Ų for all replicas and d_m values were chosen to be 11.0, 11.5, 12.0, 12.25, 12.5, 12.75, 13.0, 13.25, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5, and 17.0 Å. Both MD nuclear temperature and electronic temperature were set to 1000 K and the former was controlled by the Nośe-Hoover chain method³⁴ in the NVT ensemble. We set the MD time step to be 0.5 fs and each pair of replicas was exchanged every 40 MD steps for the former simulation and 20 MD steps for the latter simulation. Therefore, the total numbers of replica-exchange trials were 2500 and 5000 in these REUS simulations. After the simulations, we computed the PMF according to Eq. (7).

Selected DFTB-optimized structures were reoptimized using the Tao, Perdew, Staroverov, and Scuseria (TPSS)³⁴ density functional theory (DFT) with Ahlrich’s split-valence all-electron basis set,³⁵ with polarization functions on all non-hydrogen atoms. The resolution-of-identity (RI) approach³⁶ was employed as implemented in the program package TURBOMOLE.³⁷ With the exception of the iron atom, which was computed in its quintet spin state, we employed the singlet spin state in all DFT calculations.

We performed the DFTB-REUS-MD simulation for 150 ps per replica. We obtained several structures as shown in Fig. 7. After the simulation, we calculated the PMF along the reaction coordinate as shown in Fig. 8. We could obtain the PMF only in the vicinity of the FePc structure in the simulation “set A.”

In the PMF in Fig. 8, the lowest free energy is at ξ = 11.0 Å and only FePc was found as shown in Fig. 7(a) in this region. This means that FePc was the most stable structure with a reaction coordinate value between 10.5 Å and 18.0 Å. There are shoulders around ξ = 12.0 Å and 16.0 Å. In the former shoulder, almost complete FePc structures, which only failed to create the last C—N bond in the internal PN, were formed as shown in Fig. 7(b). In the latter shoulder, which is located between ξ = 15.0 Å and 17.0 Å, there were some PN trimers and these trimers had very similar structures as shown in Fig. 7(c).

Around ξ = 13.5 Å, there is a very high barrier and this barrier separates the FePc and a trimer similar to FePc. Therefore, we propose that the addition of the last PN to the oligomerized trimer is a bottleneck in the FePc formation mechanism. We also found a specific tetramer, shown in Fig. 7(b), which is very similar to FePc in this region. As we mentioned above, internal and external bonds between N and C atoms were made by the process of forming FePc. By contrast, an external bond between C and C, which should not be made by the process of forming FePc, appeared in this structure as shown in Fig. 7(d), and the mismatched bond pattern results in this structure having a reactive nitrene group that is prone to further chemical reactions. In the following, we will refer to this structure as “Structure M.” Structure M had the reaction coordinate value between ξ = 15.0 Å and 17.0 Å.

We also performed the DFTB-REUS-MD simulation for 50 ps per replica using the initial coordinates in Fig. 4, which set the four PNs and iron atoms more randomly than the initial structures in the “A set” simulations. In this simulation, we again succeeded in the formation of FePc, which had not been possible using conventional DFTB-MD simulations without REUS. We also found Structure M again, with its characteristic C—C bond and nitrene reactive center. These two structures are shown in Fig. 9. Furthermore, we found dimer and trimer structures, which had a bond between C and C like Structure M, and these dimers and trimers formed Structure M finally. We call these structures the “C—C dimer” and “C—C trimer.” The movies of forming FePc in one of the replicas (Replica 5) and Structure M in one of the replicas (Replica 2) by the “B set” REUS simulation are given in the supplementary material.

In order to study the stability of the structures, we performed DFTB optimizations of several structures, which are the initial molecular structure, PN dimer, C—C dimer, PN trimer, C—C trimer, PN tetramer, Structure M, and FePc as observed in the “B Set” DFTB-REUS-MD simulation. After eight structures were fully optimized by DFTB, we re-optimized them using the RI-TPSS-D3/def-SV(P)^35–37 method. The optimized structures from both DFTB and DFT are shown in Figs. 10 and 11 with the chemical structure files (xyz) of these optimized structures as given in the supplementary material. We also list the ΔE values relative to the energy of the FePc optimized structure at the respective level of theory in Table I for both DFTB and DFT methods.

In Table I, FePc has the lowest energy and is the most stable structure. Surprisingly, by DFTB optimization, Structure M, which has a C—C bond, is more stable than the PN tetramer. This means that this C—C bond structure is one of the metastable states in the process of forming FePc. The overall exothermicity of the reaction, starting from 4 PN molecules and one iron atom to FePc, is 260 kcal/mol at the DFTB level of theory and 308 kcal/mol at the DFT level of theory. Intermediate structures are difficult to describe at the DFT level, as the ground state is an open-shell singlet, and hence the close-shell singlet DFT energies reported here are much higher than the Mermin free energies of DFTB. Nevertheless, the DFT geometry optimizations confirm the existence of the intermediate structures, except for the PN tetramer structure, which rearranged into the FePc lower energy structure.

We further analyzed the structures, which were obtained by the “set B” DFTB-REUS-MD simulation, and found specific processes, basically following similar pathways as for “set A.” When PNs and the iron atom formed a FePc, each molecule was bound on specific pathways as shown in Fig. 12. As the first step, four PNs moved toward the iron atom and then PNs were bound with the iron atom [see Fig. 12(a)]. In the next step, carbon and nitrogen atoms in one PN, which was bound with the iron atom, created internal bonds. At the same time, this chemically altered PN, was bound with another PN and formed a “PN dimer” [see Fig. 12(b)]. In the same way, this dimer made an internal N—C bond and then formed “PN trimer” in turn by binding with the other PN [see Fig. 12(c)]. In the end, the last PN was bound with the PN trimer and a FePc was formed [see Fig. 12(d)].

We succeeded in observing the formation of FePc from one iron atom and four PNs in DFTB-REUS-MD simulations, irrespective of the initial geometries, at 1000 K. We calculated the PMF for the reaction coordinate used for the umbrella sampling and obtained interesting structures besides the thermodynamically most stable FePc, most prominently Structure M, which is a metastable structure featuring a reactive nitrene unit. Its prevalence indicates that such misbonded structures are frequently and irreversibly formed, preventing the formation of FePc and thereby causing the experimentally observed incomplete reaction with substantial byproduct in the one-pot PN method.

Moreover, by analyzing the structures and PMF, which were obtained by two sets of simulations differing only in initial geometries, we found specific pathways for the formation of a FePc from four PN and one iron atom. Based on our observations, we propose a “step-wise” pathway of forming FePc in this system and this stepwise pathway proceeds as shown in Fig. 13.

In the first step, four PNs approach the iron atom and are bound with the iron atom. In the second step, each PN connects with a neighboring PN stepwise but only three of the four PNs form a trimer. In these two steps, the reaction of forming FePc proceeds smoothly. However, the last step requires more time since the final fourth PN molecule needs to rotate and “squeeze” into the same plane as that of the PN trimer to finally form the prominent FePc structure. On the other hand, during the delicate cyclization reaction, if a bond between carbon and carbon is made, Structure M is formed. Once Structure M is formed, it is very difficult to break the strong C—C bond. Experimentally, it is known that the yield of FePc from PNs and iron or iron metal is not very high. Therefore, we also propose that Structure M is one of the final products not forming FePc.

See supplementary material for Movie 1: The movie of FePc formation from the initial coordinate by the “B set” DFTB-REUS-MD simulation in Replica 5. The frame rate of this movie was set to be 1 frame per 200 MD steps. Movie 2: The movie of Structure M formation from the initial coordinate by “B set” DFTB-REUS-MD simulation in Replica 2. The frame rate of this movie was set to be 1 frame per 200 MD steps. Also included are chemical structure files: PN_dimer_DFTB (the optimized Cartesian coordinate file of PN dimer of “B set” simulation by DFTB optimization); CeC_dimer_DFTB (the optimized Cartesian coordinate file of C-C dimer of “B set” simulation by DFTB optimization); PN_trimer_DFTB (the optimized Cartesian coordinate file of the PN trimer of “B set” simulation by DFTB optimization); CeC_trimer_DFTB (the optimized Cartesian coordinate file of the C—C trimer of “B set” simulation by DFTB optimization); PN_tetramer_DFTB: (the optimized Cartesian coordinate file of the PN tetramer of “B set” simulation by DFTB optimization); Structure_M: (the optimized Cartesian coordinate file of the C—C tetramer of “B set” simulation by DFTB optimization); FePc_DFTB: (the optimized Cartesian coordinate file of FePc of “B set” simulation by DFTB optimization); PN_dimer_DFT: (the optimized Cartesian coordinate file of the PN dimer of “B set” simulation by DFT optimization); CeC_dimer_DFT: (the optimized Cartesian coordinate file of the C—C dimer of “B set” simulation by DFT optimization); PN_trimer_DFT: (the optimized Cartesian coordinate file of the PN trimer of “B set” simulation by DFT optimization); CeC_trimer_DFT: (the optimized Cartesian coordinate file of the C—C trimer of “B set” simulation by DFT optimization); Structure_M: (the optimized Cartesian coordinate file of the C—C tetramer of “B set” simulation by DFT optimization); FePc_DFT: (the optimized Cartesian coordinate file of FePc of “B set” simulation by DFT optimization).

Some of the computations were performed on the supercomputers at the Institute for Molecular Science and at the Information Technology Center, Nagoya University. This work was supported, in part, by the Grants-in-Aid for Scientific Research (S) (No. JP16H06353), (A) (No. JP25247071), and the Grants-in-Aid for Scientific Research on Innovative Areas “Dynamical Ordering and Integrated Functions” (No. JP25102009). S. Irle was supported in part by the Laboratory Directed Research and Development (LDRD) Program of Oak Ridge National Laboratory. ORNL is managed by UT-Battelle, LLC, 651 for DOE under Contract No. DE-AC05-00OR22725.