extfieldsofstudy: Chemistry Physics Biology

s2fieldsofstudy: Chemistry Physics

text:
Methodology of Parameterization of Molecular Mechanics Force Field From Quantum Chemistry Calculations using Genetic Algorithm: A case study of methanol

In molecular dynamics (MD) simulation, force field determines the capability of an individual model in capturing physical and chemistry properties. The method for generating proper parameters of the force field form is the key component for computational research in chemistry, biochemistry, and condensed-phase physics. Our study showed that the feasibility to predict experimental condensed phase properties (i.e., density and heat of vaporization) of methanol through problem specific force field from only quantum chemistry information. To acquire the satisfying parameter sets of the force field, the genetic algorithm (GA) is the main optimization method. For electrostatic potential energy, we optimized both the electrostatic parameters of methanol using the GA method, which leads to low deviations of between the quantum mechanics (QM) calculations and the GA optimized parameters. We optimized the van der Waals (vdW) parameters both using GA and guided GA methods by calibrating interaction energy of various methanol homo-clusters, such as nonamers, undecamers, or tridecamers. Excellent agreement between the training dataset from QM calculations (i.e., MP2) and GA optimized parameters can be achieved. However, only the guided GA method, which eliminates the overestimation of interaction energy from MP2 calculations in the optimization process, provides proper vdW parameters for MD simulation to get the condensed phase properties (i.e., density and heat of vaporization) of methanol. Throughout the whole optimization process, the experimental value were not involved in the objective functions, but were only used for the purpose of justifying models (i.e., nonamers, undecamers, or tridecamers) and validating methods (i.e., GA or guided GA). Our method is general and can be extended to any type of small organic molecule systems, as well as other descriptive polarizable force field.

INTRODUCTION
Molecular mechanics (MM) modeling have been extensively employed in the research of chemistry, biochemistry, and condensed-phase physics, due to its capability of providing atomistic resolution data within relatively small to moderate computational cost. The accuracies of the MM calculation results relies on several important factors, such as sampling techniques, the underlying function forms and the parameters employed in the potential, with the latter two being the indispensable components of MM force field. Force field determines the capability of a certain model in capturing physical and chemistry properties in molecular dynamics (MD) simulations. The construction of a force field and its parameterization are non-trivial processes that require thoroughly understanding of the underlined physics and chemistry, carefully designing of parameterization protocols and appropriately applying of the optimization techniques. While several force field parameterization strategies have been widely employed and well-documented elsewhere, 1 these approaches rely more or less on both experimental and quantum mechanics (QM) calculated properties. The present study intends to extend the existing scopes of force field developing strategies, to explore the possibility to use only QM calculation approaches for force field parameterization.
In this study, we probed the procedure for parameterization of the polarizable force field in the form of Atomic Multipole Optimized Energetics for Biomolecular Applications (AMOEBA) based on QM calculations. Comparing with non-polarizable force fields, which describes the system dipole moments in an averaged and thus fixed fashion, polarizable models were designed in a way, so that they can effectively describe the capability of chemistry species to change electronic distribution under the influence of its electrostatics environment. 2 The polarizable models represented a significant improvement in the field of theoretical chemistry since last decade and had led to important scientific insights such as in the study of ionic liquids and ligand binding to biological macromolecules. 3 Knowledge on efficient and accurate parameters development of polarizable models will lead to more discoveries in biophysics, biochemistry and material science. There are several polarizable force fields available for simulating various chemical and biomolecular systems, such as the Drude polarizable force field that utilizes charges on springs, 2a, 2b, 4 the CHARMM charge equilibration (CHEQ) model that employs fluctuating charges, 5 the AMBER ff02 using partial charges with inducible dipoles, 6 and the AMOEBA force field that incorporates contributions from monopoles (charges), dipoles and quadrupoles, that are of higher order multipoles. 2c, 7 The AMOEBA polarizable force field has illustrated its capability of providing reasonable predictions of interaction energy and structural properties of small molecules and clusters in gas phase, as well as the physicochemical properties of bulk phase, such as density, self-diffusivity and static dielectric constants. 3e In order to improve the ability to describe the physical models, AMOEBA polarizable force field invokes the utilization of a vast number of variables at the atomic site. The increments of the dimension of variables lead to a series of challenges in parameter searching and thus have been a focus of the present study. To address this issue, we implement the genetic algorithm (GA) to approach the global optimization regarding parameterizations. GA is an evolutionary algorithm that mimics the process of natural selection. 8 For traditional parameters optimization, gradient-based algorithms, such as the steepest descent method, conjugate gradient method, etc. are the standard approach. 9 However, for the gradient-based method, the optimal result depends strongly on the existing knowledge of the objective function. For example, the convergence of the optimal parameters is not guaranteed when the surface of the objective function is considerably rugged or non-differentiable, or when the initial value deviates significantly from the global minimum. In parameters searching of force field, the value range of parameters exponentially expands as the complexity of the force field functional form increases. To tackle the optimization algorithm issue, statistical search heuristic (such as GA, simulated annealing, etc.) is much needed to drive the parameter space towards the optimal region for the objective function. One successful example is using the GA to determine the parameters for the reactive force field (ReaxFF) for materials system. 10 Another issue to address in this study is the selection of reference data. Force fields that have been extensively employed and validated are often based on reference data from both experiments and theoretical calculations. For example, the development of almost all the widely employed water models 2c, 11 are based on achieving a balance between the high level ab initio QM calculated dataset and condensed phase properties of bulk water measured by experimentalists. Calibration of force field parameters can be achieved by minimizing the difference of the desired properties between the MM calculations and QM calculations/experimental data. 12 For small molecules and clusters in gas phase, QM calculations are often straightforward and within reasonable computational cost. While the evaluation of the bulk phase properties from QM calculations is often computationally unfeasible, thus referring to experimental data (such as densities, dielectric constants, etc.) is necessary. However, for many molecules and especially mixtures system, experimental data are either unavailable or have significant discrepancies between different studies, due to difficulties in achieving precise measurements. Thus it is necessary to explore the possibility to select reference data of small molecules or clusters solely from theoretical calculations, i.e. ab initio QM calculations, as the training dataset for force field development to make the predictive calculation in computational modeling.

METHODS
The principal goal of this work is to establish guidance for constructing a framework of force field parameters development using primarily QM calculation data. In this initial attempt, we selected liquid methanol as the model system. Methanol represents a molecule of significant importance to organic liquids and polymers that contain bulky alkyl moieties. The hydroxyl group leads the molecule to the polar nature. The formation of the methanol clusters primarily comes from the contribution of hydrogen bond networks. 13 These ubiquitous aspects make methanol as an excellent candidate for designing and examination of force field development.
Ren et. al. suggested a general protocol of parameter development for the AMOEBA polarizable force field. 7,14 In the AMOEBA model, the electrostatic potential explicitly considers the effect of atomic multipoles (charge, dipoles, and quadrupoles) and polarization effect on induced dipole. The van der Waals (vdW) term is modeled as a buffered 14-7 potential proposed by Halgren. 15 Compared to the 12-6 potential, the buffered 14-7 potential has the advantage of optimizing structures with initial crude geometries. For hydrogen atoms, an additional parameter (the so-called reduction factor) is used to scale the position of the hydrogen atom interacting site along the corresponding covalent bond. It is meant to reflect the degree of which hydrogenic electron density displaced toward the heavy atom when covalent bonding takes place. 3e The full potential energy function of the AMOEBA force field has been described in full detail elsewhere. 2c, 14 The general energy expression of the AMOEBA force field is written in equation (1), where parameterizations of electrostatic potential energy E ESP and vdW interaction energy E vdW terms are studied in this work. Bonded interaction, including bond, angle, dihedral and torsional parameters, were kept consistent with that in the original AMOEBA force field (amoeba09.prm). 7 In Ren et. al.'s protocol, the parameters of multipoles and torsional bonded terms are optimized using QM results. However, besides of taking QM calculations as the reference, the protocol adapts data from existing database of experimental values, which are not entirely available for every desired system. We followed their approaches with the context of parameterizing standard AMOEBA force field for atomic multipoles, polarizabilities and vdW interaction by GA methods using solely results from QM calculations, which are available for the methanol system within given computational resources.
All the QM calculations were performed with Gaussian 09 16 package of electronic structure programs. The second order Møller-Plesset perturbation theory 17 (MP2) was employed, with the basis set superposition error (BSSE) correction. 18 This level of theory and basis set to predict energetics and structural properties has been widely verified by previous studies. 19 The electronic density calculation of methanol monomer is initially used for calibrating the magnitude of the multipole moments at the atomic positions.
Combining with the Distributed Multipole Analysis (GDMA 2.2) tool developed by Anthony Stone, 20 the Tinker 2c, 21 package is used to optimize the permanent atomic multipoles by fitting the electrostatic potential measured on the Connolly surface of the methanol monomer. Table S1 in the Supporting Information is showing the electrostatic potential results for methanol monomer using different basis sets.
In the parameterization process, the training dataset for the optimization should contain a good sampling of the possible structural configurations and their respective energies. For instance, to predict the polarizable effect between molecules, methanol dimers with a continuous range of electrostatic potential energies ( E ESP ) should be included. In AMOEBA, to compute clusters electrostatic potential energy ( E ESP ) including the polarizable effect between methanol molecules, 44 independent parameters are needed. Those parameters are monopole (q), dipole (µ x , µ y , µ z ), quadrupole-a traceless and systemic matrix-(Q xx , Q yx , Q yy , Q zx , Q zy ), atomic polarizability (α) and Thole's description of damping factor (a) for four types of atoms (O, H(-O), C and H(-C)) in methanol (The detailed description of the parameter form can be found in somewhere else 2c, 14 ). We optimized all 44 independent parameters listed above, using an extensive training dataset taking electronic density calculations for 4943 methanol dimers. The methanol dimer configurations were optimized via Guassian09 program at the MP2 level using 6-311G(d, p) basis set. The electrostatic potential energy ( E ESP ) of 4943 dimers was computed using the same level of theory and basis set. The detailed description of optimization for electrostatic parameters is in the Supporting Information Section S2.
To ensure an adequate representation of the various possible coordination environments and cluster sizes as well as the energy landscape, we sampled three training datasets containing different methanol homoclusters, which consisting of 9, 11, and 13 molecules referred as nonamer, undecamer and tridecamer in the descriptions below, respectively, to calibrate the vdW interaction energy ( E vdW ). For methanol, there are vdW potential depth (ε) and minimum energy distance (R * ) for four types of atoms, and additional reduction factor (λ) for two types of hydrogen atoms, which constitutes 10 independent parameters to optimize. We optimized these 10 independent parameters using interaction energy ( ΔE ) results of the sampled 502 nonamers, 157 undecamers and 94 tridecamers at the MP2/6-31G (d, p), respectively. For instance, the nonamer was sampled through randomly placing 9 C atoms in a 10 Å × 10 Å × 10 Å computational supercell, where one C atom placed in the center and surrounded by other 8 C atoms with relative distance at least of 3.25 Å, the minimum distance of carbon bond networks in methanol, 22 to each other. Then, methanol with random orientations will be assigned to the position where the C atoms sit. The methanol nonamer configuration is optimized in Gaussian09 using the eigenvalue-following algorithm 23 with constraints of fixed C atoms positions. The relaxed nonamer configurations and their interaction energies ( ΔE ) between the center methanol and surrounding methanols are then employed in the training dataset for optimizing the vdW parameters. This procedure was applied to undecamers and tirdecamers as well for getting the relaxed configuration and the interaction energy, respectively. Figure 1 (a), (b) and (c) are showing example configurations of nonamer, undecamer and tridecamer. In the AMOEBA force field, except for the electrostatic potential energy ( E ESP ), vdW interaction energy ( E vdw ) is the only term consisting of interaction energy between molecules. The vdW interaction energy ( E vdw ) can be calibrated from the interaction energy ( ΔE ) by separating apart molecules, with the assumption that the deformation energy of each molecule upon binding is reasonably small. 7 Through sampling configurations widely and calculating the coordinated interaction energies, we consider the vdW interaction by taking the effect of electrostatic interaction inclusively. The interaction energy ( ΔE ) is calculated as the total energy of the cluster subtracts the total energy of the center methanol molecule and that of the surrounding methanol clusters, as shown in equation (2). The BSSE correction is applied for getting the ΔE of every configuration. The MP2/6-31G (d, p) computed interaction energies of the 502 nonamer, 157 undecamer and 94 tridecamer configurations exhibit a Gaussian distribution as shown in figure 1 (d), (e) and (f), which signify the sampling of those configurations are adequate and in the equilibrated states.
Using the training dataset described above, we optimized the non-bonded parameters in AMOEBA force field by employing genetic algorithm; the procedure is outlined in Figure 2. Taking the optimization of vdW parameters as an illustration, we begin the optimization process by generating a population of N p = 120 parameter sets randomly, such that their values lie within physically allowable limits (parameter search ranges are listed in Supporting Information Table S2). Each set of parameters in this population is called a member. For each member i, we compute the interaction energies for all structures in the training dataset using MD simulation package Tinker 2c, 21 and evaluate the objective function Δ i given by equation (3), where ΔE j MP 2 and ΔE j AMOEBA i are the MP2 calculated and the AMOEBA calcuated interaction energies for the structure j in the training dataset, respectively. Since the electrostatic potential energy ( E j ESP ) is computed for the structure j to represent the atomic electrostatic potential precisely, the only uncertain term in the interaction energy is the vdW interaction energy ( E j vdW ), represented by a member i, a set of 10 independent parameters. Here, members are then ranked in an ascending order of Δ i . After the ranking, non-linear roulette wheel selection 24 will be performed to select the top 60% members, i.e., the ones with lowest values of Δ i , which are then subjected to genetic operations: crossover with crossover-rate 3% and mutation. These mutations introduce sufficient diversity into the population, and the non-linear selection scheme helps to avoid premature convergence of the GA run. After the genetic operations, both the old and the new members are ranked by their Δs. The best N p parameter sets (members) were then chosen to constitute the next generation. Such an optimization routine ensures that only satisfactory parameter sets survive after each generation; upon repeating this workflow for sufficient generations and sampling viable regions in the parameters space, we perform three separate GA runs starting with different random populations. From each of the converged GA run, we choose the final parameter set corresponding the lowest Δ.
To overcome one of the serious shortcomings of MP2 theory, a noticeable overestimation of the interaction energy, 25 which plays a major role in stability and formation of condensed phase molecular structure, we improved the form of objective function Δ in the GA program. We introduced a penalty parameter δ in the GA optimization procedure as the unexplained discrepancy between the AMOEBA results and MP2 calculations, as shown in equation (4), which obeys the design concept of guided genetic algorithm (guided GA). 26 This penalty parameter δ never participates as the form of AMOEBA force field in MD simulations, adds up to 11 independent parameters as a subset of i th member in the GA optimization process for vdW parameters.

RESULTS
The vdW parameters are critical in the force field for getting the correct condensed phase structures, such as right intermolecular distance, through molecular dynamic simulations.
In this section, we present the optimization results of the vdW parameters from both GA and guided GA program by calibrating interaction energy between molecules. We follow with showing binding curves for methanol dimers using both GA and guided GA optimized parameters. Finally, we show the MD simulation results (i.e., density and heat of vaporization) of a condensed phase methanol system from both GA and guided GA optimized parameters using AMOEBA force field.

VDW parameters of nonamers, undecamers, and tridecamers from GA and Guided GA
The GA and guided GA programs were applied to optimize the vdW parameters using the interaction energies ( ΔE ) calculated at MP2/6-31G (d, p) with BSSE correction for 502, 157 and 94 methanol nonamers, undecamers and tridecamers, respectively. The excellent correlations (R = 0.970, 0.976, and 0.968) of the ΔE between MP2 calculations and AMOEBA with GA optimized parameters for nonamers, undecamers, and tridecamers, respectively, are shown in figure 3 (a), (b) and (c). Meanwhile with apparent overestimation of the MP2 calculation, figure 3 (d), (e) and (f) are also showing excellent correlations (R = 0.956, 0.965 and 0.946) of the ΔE between MP2 calculations and AMOEBA with guided GA optimized parameters for nonamers, undecamers, and tridecamers, respectively. There are in total six sets of vdW parameters from nonamers, undecamers, and tridecamers all using GA and guided GA optimization method, respectively.

VDW parameters from Dimers
To validate the transferability of the GA and guided GA optimized AMOEBA parameters from nonamers, undecamers, and tridecamers for gas-phased methanol molecules, we use the six sets of parameters to calibrate the binding energy curves of methanol dimers. We fully optimized 141 methanol dimer configurations at MP2/6-31G (d, p) with BSSE correction with the constraint of pre-defined distance between carbon atoms in the two methanol molecules. Figure 4 is showing the binding energy of methanol dimers as a function of the distance between two carbon atoms.
We obtained the binding energy curves for the 141 configurations of methanol dimer using the three sets of vdW parameters from nonamers, undecamers, and tridecamers using GA optimization method, respectively. Figure 5 (a), (b) and (c) are showing the good agreement of the binding energy between the MP2 calculation and the GA optimized parameters from nonamers, undecamers, and tridecamers, respectively. We also obtained the binding energy curves for those dimers using the other three sets of vdW parameters using guided GA optimization method from nonamers, undecamers, and tridecamers. Figure 5 (d), (e) and (f) are showing the binding energy curves of methanol dimers from the MP2 calculations and from the guided GA optimized vdW parameters by nonamers, undecamers, and tridecamers, respectively. From the lower row of figure 5, we can see discrepancies of binding energy between the gas-phased MP2 calculations and the guided GA optimized parameters from clusters, where the overestimation of MP2 calculations is conspicuous near the lowest energy curve. The shifting down of the binding energy curves from the guided GA optimized parameters implies the interaction energy of the same level of QM theory (i.e., MP2/6-31G(d, p) + BSSE) calculations from gas-phased molecules leads to the overestimation of real potential energy surface (PES). This discrepancy of interaction energy from MP2 calculations is commonly seen as the incompleteness of the basis sets for the particular level of MP2 theory and the insufficiency sampling of configurations in the parameterization of force field. 19,25

VDW parameters for Condensed Phase Methanol
To verify the ability of our parameters of predicting the condensed phase properties, sets of MD simulations using the developed parameters of AMOEBA force field were performed. Here, we applied the six sets of optimized parameters to perform MD simulations to the condensed phase system consisting of 344 methanol molecules. For getting the density (ρ) and heat of vaporization (H vap ), we used Tinker package (omm_dynamics program --GPU version of Tinker dynamics program 27 ) to run MD simulation. NPT ensemble simulations were performed at T = 298.15 K and P = 40 atm with an integration time step of 1 fs. With three different random seeds for initiating velocity on the system, the averaged value of density and heat of vaporization were obtained. Table 1 shows the result of density and heat of vaporization from Tinker MD simulations with the GA optimized AMOEBA force field parameters. Among all the results, the simulated density from tridecamers is the best, while the winner of the heat of vaporization is from undecamers. However, from Table 1, it is clear that the value of density and heat of vaporization are nowhere close to the experimental value.  Table 2 shows the results of density and heat of vaporization from Tinker MD simulations with the guided GA optimized AMOEBA force field parameters. The simulated density and heat of vaporization from tridecamers are much closer to the experimental values. We can see that using the guided GA, which introduces a penalty parameter δ, eliminates the overestimation of interaction energy from MP2 calculation. The guided GA optimized parameters describe better condensed phase properties (i.e., density and heat of vaporization) than the GA optimized parameters in MD simulation. Among all the results, the simulated density and heat of vaporizations from more number of methanol molecule clusters are better. MP2 calculations of interaction energy of tridecamers and undecamers are more descriptive than nonamers for optimizing vdW parameters, due to more detailed information can be calibrated as more molecules in the solvation shell of methanol involved.

CONCLUSIONS
Our study showed that the feasibility to predict experimental condensed phase properties (i.e., density and heat of vaporization) of methanol through problem specific (here, we use AMOEBA) force field from only quantum chemistry information. To acquire the satisfying parameter sets, the genetic algorithm (GA) is the main optimization method. For electrostatic potential energy ( E ESP ), we optimized the both the multipoles and polarizability, Thole damping factors of methanol using the GA method, which leads to lower deviations of E ESP between QM calculations and the GA optimized parameters than that from amoeba09.prm. We optimized the vdW parameters both using GA and guided GA methods by calibrating interaction energy ( ΔE ) of various methanol homo-clusters, such as nonamers, undecamers, or tridecamers. Excellent agreement between the training dataset from MP2 calculations and GA optimized parameters can be achieved. However, only the guided GA method, which eliminates the overestimation of interaction energy from MP2 calculations in the optimization process, provides proper vdW parameters for MD simulation to get condensed phase properties (i.e., density and heat of vaporization) of methanol. Throughout the whole optimization process, the experimental value were not involved in the objective functions, but were only used for the purpose of justifying models (i.e., nonamers, undecamers, or tridecamers) and validating methods (GA or guided GA). We conclude that the main difficulty of parameterizing the force field for liquid phase methanol from solely QM calculations could be coming from the reliability of the training dataset (i.e., MP2 calculations), as it is well known that different basis sets of MP2 could lead to different accuracy of QM calculations, the conspicuous overestimation of interaction energy from MP2 calculations. Some of the other higher level of QM calculations (e.g., MP3/4, CCSD(T), etc.) can only deal with up to certain amount of atoms/electrons. 29 It is ambiguous to study the property of condensed phase properties from a few molecules. In fact, we tried to use the interaction energy ( ΔE ) calculated for methanol tridecamers from another type of QM calculation, Symmetry-Adapted Perturbation Theory (SAPT), as the training dataset using GA optimization, which yields a lousy performance of the condensed phase properties. However, our method was able to utilize the QM calculations (i.e., MP2 calculations) from the limited calculation resource to tune the force field parameters for reasonable agreement between the MD simulated values and experimental data. The method is general and can be extended to any type of small organic molecule systems, as well as other descriptive polarizable force field.


Transmission and Reflection Terahertz Spectroscopy of Insensitive Melt-Cast High-Explosive Materials

Currently, artillery shells and grenades that are introduced into the market are based on melt-castable insensitive high explosives (IHEs), which do not explode while they run a risk of impact, heat or shrapnel. Particles of explosives (such as hexogen, nitroguanidine and nitrotriazolone) are suspended in different proportions in a matrix of 2.4-dinitroanisole. In this paper, we investigated samples of commonly used IHEs: PAX-41, IMX-104 and IMX-101, whose internal structures were determined by a scanning electron microscope. Terahertz time domain spectroscopy was applied in both transmission and reflection configurations. At first, the complex refraction indices of four pure constituents creating IHEs were determined and became the basis of further calculations. Next, the experimentally determined transmission and reflection spectra of IHEs and pure constituents were compared with theoretical consider-ations. The influence of the grain size of constituent material and scattering on the reflection spectra was analysed, and good agreement between the experimental and theoretical data was achieved.

1 Introduction materials manufactured on the basis of TNT may easily detonate because of impact with shrapnel, falling or temperature escalation. Bearing in mind numerous cases of accidental human victims over the last several years, mixtures based on TNT are replaced by insensitive high explosives with 2.4-dinitroanisole (DNAN) as the melting ingredient. Recently, the US Army introduced a family of melt-castable DNAN-based explosives containing mainly mixtures of RDX, 1-nitroguanidine (NQ) and 3-nitro-1,2,4-triazole-5-one (NTO) [1][2][3][4]. Trace amounts of N-methyl-p-nitroaniline (MNA) are included in some formulations to aid in the processing. IMX-101 (mixture of DNAN, NTO and NQ) replaces TNT as the main filler explosive in artillery projectiles. IMX-104 (mixture of DNAN, NTO and RDX) is qualified as a replacement for Composition B in mortars. PAX-41 containing DNAN and RDX is tested for spider grenades [1][2][3][4].
More and more growing threats caused by the use of explosive materials are forcing constant development of novel and modern methods for their detection [5]. One of the driving forces behind the development of the terahertz technique (the spectral range from 0.1 to 10 THz) is the fact that many commonly used explosive materials (such as RDX, PETN, HMX, TNT and NQ) have unique spectral features in this range of radiation [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Moreover, non-polar dielectric materials (such as clothes, paper and plastics) are relatively transparent in the range up to about 3 THz [6,9], which raises hope of efficient identification of explosive materials remaining under cover or being inside packaging [6,12]. Spectral examinations are mostly carried out by means of time domain spectroscopy (TDS) [6-11, 13, 17-19], Fourier spectrometry [10] and the use of a tunable source and broadband detector [12]. Most of the mentioned explosive materials have been thoroughly investigated in transmission [6-8, 13, 16] and in reflection configurations [7,13]. Experiments have been carried out in different temperatures [14]; also, solid-state density functional theory calculations have been performed [15]. Methods for identifying explosive materials by means of their spectra are constantly being developed [16][17][18][19]. RDX is the most widely investigated material due to the fact that it is a constituent of many commonly used plastic materials (such as C-4 and Semtex) and meltcastable materials (like Composition B). RDX has a very clear absorption peak over relatively low frequencies-about 0.82 THz. This peak is crucial for its identification and is accompanied by a few other less clear peaks over higher frequencies [6-8, 10, 11, 13]. NQ has been quite well described in the literature as well-it has strong absorption peaks at 1.43, 2.38 and 2.76 THz [6]. In contrast, DNAN has only been briefly characterized [20]. Terahertz investigations of NTO and IHEs were not found at all in the literature.
First, it should be noted that sample preparation is a very important factor in obtaining reliable spectra in the terahertz range. The commonly used transmission configuration [6-8, 13, 16] gives clearly visible absorption peaks that may later be used as the basis for the identification of explosive materials. With regard to quite big attenuation (∼100-300 cm −1 ) [6], practical applications are limited to rare situations such as measurements of the thin layer of the powder placed in the envelope. For this reason, all samples prepared to transmission measurements are in the form of pressed pellets consisting of material with good transmission in the terahertz range (such as polyethylene) with suspended explosive material particles with a concentration of 5-20 %, which allows for proper determination of absorption characteristics [6] but ensures sufficient transparency. In real applications, only reflection configuration can be applied because terahertz radiation cannot transmit through targets thicker than a few millimetres. Samples for reflection measurements are mostly manufactured in the form of pressed pellets but consisting of pure explosive material powder, sometimes with a few per cent addition of desynthesizer [7]. The interaction of the terahertz wave with the surface of the sample depends on its morphology. Specular reflection takes place when the sample has a flat and smooth surface. If the particle size of the solid material is comparable to the terahertz wavelength, the scattering may cause a significant change in the obtained spectrum [21]. First, we manufactured pressed pellets of considered explosive materials. To approach the most realistic conditions, we also prepared melt-castable samples, where the solubility of materials in DNAN and grain sizes of particles were analysed.
Here, we investigated three prospective IHEs, PAX-41, IMX-104 and IMX-101, by means of a terahertz TDS system both in transmission and reflection configurations. Samples were prepared as pressed pellets and melt-castable discs. Their surface morphology and grain sizes were determined by means of a scanning electron microscope (SEM). For terahertz measurements, the TDS method was used due to its sensitivity and high measurement dynamics. First, the complex transmission characteristics of four materials (RDX, NQ, NTO and DNAN) were determined to provide the exact data for further calculations. Next, the transmission characteristics of pellets manufactured from IHEs and the reflection characteristics of pellets from pure materials and discs made from IHEs were measured. The results of IHE measurements, both in transmission and reflection configurations, maintain good agreement with theoretical calculations. Images from SEM investigations help to draw conclusions about the phenomena taking place during the reflection.

Sample Preparation and Characterization
Although the literature concerning IHEs is quite impressive, the details of their manufacturing process are not available due to the sensitive character of such data and patent claims. Because of this, samples were prepared based on the literature [1][2][3] as well as on our experience. Three types of samples were prepared: pressed pellets for transmission (denoted by the suffix p) and reflection (denoted by the suffix P) measurements and melt-castable discs (denoted by the suffix C) for reflection measurements ( Table 1, Fig. 1).
For transmission configuration, the explosive materials were ground using a mortar and pestle to reduce the particle size below 50 μm, which suppresses scattering losses [6,20]. Due to the high attenuation of pure explosives, we mixed them with fine polyethylene (PE) powder (used as a matrix), which is highly transparent in the terahertz range [6]. The mixing ratio was 5 % weight of a material to total weight of mixed samples. Afterwards, the resultant fine powder was pressed into a pellet directly using a hydraulic press. As a result, we obtained the 13-mm-diameter pellets listed in Table 1 (numbers 2-8), which consisted of 20 mg of explosive material and 360 mg of PE (Fig. 1a). As a reference, a 360-mg pellet pressed from pure polyethylene was prepared (number 1). Although the powders were thoroughly mixed, the low concentration of the constituent materials in the sample (e.g. 1.25 % RDX in IMX-104) can result in fluctuations in their distribution, which may later lead to disturbances of the obtained terahertz spectra. For reflection investigation, fine powders of pure explosives and IHE materials (numbers 9-15) were pressed into pellets with a weight of about 5 g, thickness of 5 mm and diameter of 30 mm (Fig. 1b). Macroscopically, surfaces of the samples were flat and smooth. However, SEM investigations (by Quanta 3D FEG) revealed that the surfaces were slightly rough and consisted of grains of crystals with dimensions below 100 μm and flat surfaces (Fig. 2a-c). In all reflection measurements, a gold mirror was used as a reference.
The melt-castable samples were prepared as 30-40-mm-diameter discs with one flat and smooth surface, which suppresses surface-dependent scatterings and allows the interaction of the radiation with the material to be focused upon (Fig. 1c). The powders of RDX, NQ and NTO (synthesized at the Military University of Technology according to the original procedures) were sieved by using a vibratory sieve shaker AS 200 control (Retsch), in order to obtain the desired grain diameters. Cast samples were prepared by using laboratory-scale meltpour apparatus equipped with a thermocycler and heated pot. DNAN was heated in the pot to a temperature of 110°C until it was totally melted. Next, a mechanical stirrer was started up and solid ingredients were added in small portions. The mixture was stirred for 15 min until the composition was homogenous. The casting was realized by tilting the melting pot, and the composition was poured into a rubber mould placed on glass. After solidification, the cast disc was recovered and tested using the methods described below. The mixture of DNAN and RDX was studied more deeply to enable understanding of solubility-related phenomena. Since the solubility of RDX in DNAN is about 10 % (at 100°C), we manufactured the D90R10 sample (number 16), which consisted of 90 % of DNAN and 10 % of RDX. Next, we prepared three samples with various mean grain sizes of RDX particles (below 100, 200-250 and 300-350 μm) to study the influence of solubility and scattering. They are called PAX-41_100, PAX-41_200 and PAX-41_300, respectively (numbers [17][18][19]. The sample of IMX-104 (IMX-101) consisted of mixed powders of RDX and NTO (NQ and NTO) with a particle   Both the internal structure of the cast disc (revealed from small broken-off pieces) and its surface were imaged using an SEM microscope. The surfaces of the manufactured samples were flat and macroscopically smooth ( Fig. 2d-f, k) due to the fact that smooth and flat glass was used in pouring compositions into rubber moulds. Small scratches, holes and micro-cracks that were caused by the solidification process had dimensions up to 10 μm and are not significant from the point of view of terahertz radiation.
The D90R10 sample should consist of totally molten RDX. But during cooling down, a part of RDX is recrystallized to the form of crystalline grains smaller than 20 μm. The exact efficiency of recrystallization and the size of those grains are not known to the authors. Assuming that 50 % of RDX will transform into a crystalline form, the D90R10 sample will consist of 95 % of solid solution of RDX in DNAN and 5 % of small grains. Figure 2l shows the homogenous structure of the D90R10 sample. In samples from the PAX-41 group, small RDX crystals are totally dissolved while bigger crystals have a reduced size. After cooling down and recrystallization, the samples consist of around 63 % of solid solution of RDX in DNAN and also of small grains of RDX and partially melted original grains. For the PAX-41_100 sample, the characteristic polyhedral RDX crystals (with dimensions of 10-80 μm) are clearly seen (Fig. 2g). For the PAX-41_300 sample (Fig. 2j), a large central crystallite (around 250 μm) is shown together with a number of smaller crystallites with sizes of 20-100 μm.
The IMX-101_C sample consists of around 37 % of solid solution of RDX in DNAN, less than 20 % of RDX crystals of different sizes, but all smaller than 150 μm, and around 40 % of oblong cuboid NTO crystals of diverse sizes up to 150 μm. Figure 2h shows a cross section of the sample, where the surface of the sample can be clearly seen and cuboid NTO crystallites located on the surface are visible together with holes after RDX crystallites. The IMX-101_C sample consists of dominant needle-shaped grains of NQ and cuboid NTO crystallites suspended in a solid solution of DNAN (Fig. 2i). The sizes of the NQ and NTO crystals are below 150 μm.

Experimental and Data Analysis
In the range of terahertz radiation, time domain spectroscopy is the most commonly applied method and has already been thoroughly described [6-11, 13, 17-19]. Here, we carried out an experiment using a TPS 3000 (TeraView) spectrometer, which is a commercially available TDS set-up ready to operate in transmission and reflection configurations. In transmission, the pellets were located in the focus of the optical system. The reflection set-up had an incident angle of α = 22° (Fig. 3) and was an extension of the transmission arrangement created by adding flat and off-axis parabolic mirrors with 50-mm focal lengths. In reflection, the sample was located on a stage and moved manually in X-Y directions to measure 20 different points on the sample. The surface of all samples and the gold mirror were positioned exactly in the focus of the optical system with an error of about 10 μm, which was controlled by the micrometer screw (in the Z direction). This condition is crucial for appropriate determination of the reflection spectrum; however, a misplacement phase error correction was also developed [22]. In the transmission (reflection) set-up, we also determined beam diameters using the knife-edge method [23]; these are equal to about 1.9, 0.9 and 0.5 mm (1, 0.7 and 0.4 mm) at frequencies of 0.5, 1 and 2 THz, respectively.
The TDS system used a pulse lasting about 0.4 ps (which is full width at half-maximum value). This signal has a 0.1-3-THz useful spectral range (Fig. 3). Each 40-ps-long waveform was measured with a resolution of 2048 points and averaged 3000 times. Then, Blackman-Harris third-order filtering was applied to apodize the time signal, which smoothed the spectrum. Therefore, 3 min was enough for a single measurement. In Fig. 3, the power spectrum of the pulses and the noise levels are shown for transmission and reflection configurations. Such spectra are available when the influence of water vapour is minimized by purging the chamber with dry air. To assess system performance, the signal-to-noise ratio and dynamic range were determined in the frequency domain according to the algorithm presented in [24] (Fig. 3).
The main advantage of the transmission TDS technique is that one can easily determine the complex refractive index (N = n + i · κ) of a measured material [6]. The power absorption coefficient (α) can be calculated as where c is the speed of light, ν is the frequency and κ is the extinction coefficient. The presented method requires precise determination of the thicknesses of the reference (PE) and the sample (with explosive) pellets.
Assuming that the s-polarized beam incidents from the air on the sample, its complex amplitude reflection coefficient (r) can be calculated from the Fresnel formula as where r is the amplitude reflection coefficient and Δϕ is the reflectance phase. The reflectance R is related to the amplitude reflection coefficient by R = r 2 . The reflectance spectrum (R) mimics the refractive index (n) variations, while the reflectance phase (Δϕ) represents the absorption features similarly to the extinction coefficient (κ) in transmission [7]. For unary pellets with small-size particles (samples 9-12), we used the N index directly taken from transmission measurements (Fig. 4). The interaction between the terahertz waves and the multi-component samples (26)(27), which are a dielectric mixture, can be quantitatively described using the effective medium theory (EMT) [25,26]. Although EMT assumes that the particle size is much smaller than the wavelength of the electromagnetic radiation, we proved that the experimental results are also consistent with the theory for the particles equal to Fig. 3 The power spectrum of the pulses and the noise levels in transmission and reflection configurations. SNR and DR denote the signal-to-noise ratio and the dynamic range, respectively. The inset shows a schematic diagram of the reflection configuration or slightly smaller than the wavelength. EMT considers that the wave is averaged over the volume of the sample and experiences an effective refractive index N eff . Due to the variety of sizes, shapes and orientations of RDX, NTO and NQ particles diluted in the DNAN (which is a uniform matrix), we calculated N eff as the average refractive index of the constituent materials (Figs. 5, 6, 7, 8 and 9), taking into account their volumetric composition. The correlation coefficient (C) was used as a parameter, which quantitatively describes the agreement between experimental and theoretical data and is included in Figs. 5, 6, 7, 8, 9 and 10.

Transmission Results and Discussion
Complex spectra of the RDX_p, NQ_p, NTO_p and DNAN_p samples (Fig. 4) reveal the absorption peaks which are connected with slopes of the refractive index characteristics. In the case of RDX and NQ, the positions of peaks agree well with the literature [6-8, 10, 11, 13]. Some discrepancies can be attributed to the different manufacturer of materials and fluctuations in their density inside the sample. NTO shows two distinctive absorption peaks at 1.54 and 2.75 THz, while DNAN has a broad and flat peak at 2.14 THz. Figure 5 plots the spectra of the PAX-41_p, IMX-104_p and IMX-101_p samples directly measured in the TDS set-up and calculated based on the spectra of their constituents (Fig. 4) and their concentration ( Table 1). The good agreement between the measurement and the calculation proved that the features appearing in the absorption and the refractive index characteristics can be understood from the sum of the constituents. Slightly shifted features of RDX (at 0.82 and 1.58 THz) and DNAN (at 2.2 THz) are observed for PAX-41, because it consisted of approximately equal amounts of constituent materials. In the case of IMX-104, the feature at 0.82 THz can be attributed to RDX, the one at 1.55 THz to RDX or NTO and the one at 2.12 THz to DNAN. For IMX-101, high NQ peaks dominate over DNAN and NTO. Some observed discrepancies were due to the local fluctuations in the concentration and thickness of material both in the case of samples with pure explosives (Fig. 4) and IHEs (Fig. 5). Figure 6 shows the complex reflection spectra (the reflectance R and the phase Δϕ) of 30-mmdiameter pellets: RDX_P, NQ_P, NTO_P and DNAN_P. The spectra were directly measured in the TDS set-up and calculated for comparison based on Eq. 2 and the measured absorption spectra (Fig. 4). As expected, the reflectance characteristic mimics the refractive index variations, while the reflectance phase is similar to the absorption characteristic. The RDX characteristic has an amplitude decrease (slope) of about 7 % centred at 0.87 THz. For NQ, we can see three distinctive features in the reflectance phase at 1.46, 2.52 and 2.98 THz that correspond to the amplitude slopes of 10, 20 and 7 % in the reflectance characteristic, respectively. The NTO spectrum contains a feature whose amplitude slope is about 1.5 % at 1.5 THz, while DNAN has a 2 % decrease centred at 2.2 THz.

Reflection Results and Discussion
The measured reflectance and phase characteristics are shifted towards higher frequencies in comparison with the calculated spectra. The phase spectrum of RDX (Fig. 6c) has a Fig. 6 The reflectance spectra of the NQ_P and RDX_P (a, c) and NTO_P and DNAN_P (b, d) samples The reflectance of the 30-mm-diameter pellets of IHE samples (numbers 13-15) depended on the point of measurement on the sample. Figure 7 presents only the best characteristics obtained for each sample. The number of points with the appropriate spectrum was equal to about 90, 40 and 50 % for PAX_41_P, IMX-104_P and IMX-101_P, respectively. For PAX-41_P (Fig. 7a), the RDX features at 0.82, 1.38 and 1.56 THz are clearly visible for two different measurement points. Scattering of the terahertz beam at the RDX particles is noticeable at higher frequencies, and the feature of DNAN at 2.15 THz is ambiguous. Probably, due to local fluctuation of constituent distribution, in the case of the IMX-104_P Fig. 8 The reflection spectrum of the D90R10 sample Fig. 9 The reflection spectra of PAX-41-calculated and measured at two points for samples with three grain sizes (100, 200 and 300 μm diameter) sample, which had a lower concentration of RDX (25 %) in comparison to PAX-41_P (40 %), the feature of RDX at 0.82 THz has a higher reflectance slope with an amplitude equal to about 3 % (Fig. 7b) in comparison to PAX-41_P (2.5 %). The features of NTO at 1.58 THz and RDX at 1.93 THz are not always present. Three NQ-related features in the reflectance spectrum of IMX-101_P are clearly visible at 1.45, 2.45 and 2.94 THz (Fig. 7c). The calculated spectra (black curves) agree well with the experiment up to about 1.5-2 THz. At higher frequencies, both the local fluctuation in constituent distribution and the scattering modified the experimental characteristics.
The exemplary reflection spectrum of the melt-castable D90R10 sample (Fig. 8) revealed a 1.5 % RDX-origin slope in reflectance at 0.81 THz conjugated with the peak of the reflectance phase (Δϕ). The feature at 2.25 THz can be attributed to both DNAN and RDX. A similar result could be obtained at almost every point of the sample, because it was uniform, since RDX was totally dissolved in DNAN. The proposed effective medium theory is consistent with the experiment. Since this sample consisted of 95 % of homogenous solid solution of RDX in DNAN and 5 % of small (below 10 μm) recrystallized RDX grains, its surface and volume (Fig. 2k, l) were uniform for the terahertz radiation with a beam diameter of about 0.5-1 mm and wavelength in the range 100-1000 μm. Some discrepancies at higher frequencies are probably the result of the scattering phenomenon.
Surfaces of the cast samples from the PAX-41 family with different grain sizes (100, 200 and 300 μm) were smooth and uniform. They consisted mainly of the homogenous solid solution of RDX in DNAN (Fig. 2d). However, the reflectance depended on the measurement point on the sample-in some of them, the RDX features at 0.82 THz (both the slope in the reflection curve and the peak in the phase characteristic) were deformed and in a few of them not visible at all (Fig. 9). Although the concentration of the RDX in DNAN was constant and Fig. 10 The reflection spectra of IMX-104_C (a) and IMX-104_C (b)-calculated and measured at two points equal to 40 %, the number of points without the appropriate spectrum were equal to about 10 % for PAX-41_100 and 20 % for both PAX-41_200 and PAX-41_300.
For PAX_100, smaller particles of RDX were probably more uniformly packed in the sample and the scattering was also lower. Figure 9 presents reflectance curves of the PAX-41 family measured at two selected points with the highest slope of RDX at 0.82 THz. For 100 and 200 μm grain sizes (Fig. 10a, b), the amplitude of the RDX slope at 0.82 THz is two times bigger than the EMT-based calculations, which may be attributed to the local concentration of RDX grains. For PAX-41_300 with higher grain size, the scattering increased, which resulted in a smaller amplitude of the slope at 0.82 THz, while the average reflectance remained constant and equal to about 0.09 (Fig. 9). The feature of DNAN at about 2.1-2.2 THz is not visible. At higher frequencies, the scattering suppresses the expected features for all cases.
The surface of the melt-castable disc, as in the PAX-41 case, was smooth and uniform and consisted mainly of the solid solution of RDX in DNAN (Fig. 2e). Moreover, the cross section of the IMX-104_C (Fig. 2h) presents a solid solution of RDX in DNAN and less than 150-μmsize particles of RDX (polyhedral) and NTO (rectangular). Since the IMX-104_C sample consists of 40 % of NTO and 25 % of RDX, the feature of RDX at 0.82 THz has a smaller reflectance slope than PAX-41 with an amplitude equal to about 1.2 %. However, features of NTO at 1.56 THz and DNAN at about 2.23 THz can be seen (Fig. 10a). Although the visibility of the features is rather small, the number of points with the appropriate spectrum is equal to about 90 %. The calculated spectrum agrees well with the experiment, but the amplitude of the RDX slope at 0.82 THz is smaller than expected, while the NTO-related amplitude at 1.56 THz is bigger, which probably resulted from fluctuation in the concentration of constituent grains.
As expected, the reflection spectrum of the IMX-101_C sample is dominated by NQ features, which are visible at 1.42, 2.37 and 2.79 THz (Fig. 10b). Due to the high scattering at higher frequencies, the feature at 2.79 THz is not always present. However, the melting process broke up some of the larger needle-shape NQ grains into smaller ones (Fig. 2i) and we think that the surface layer of the sample consisted mainly of a mixture of melted DNAN with fine NQ grains, which justifies the very good visibility of the NQ features even at higher frequencies. At 1.42 and 2.37 THz, the reflectance slope has an amplitude equal to about 1 and 5 %, respectively. The features of NTO at 1.54 THz and DNAN at about 2.1 THz are not visible. Probably due to clustering of the NQ grains, the number of points with an appropriate spectrum is about 40 %. The EMT-based calculated spectrum is in line with the experiment.

Summary and Conclusions
The examination of explosive materials by means of terahertz transmission and reflection spectroscopy is very promising. In the transmission configuration, the absorption peaks are connected with the slopes of the refractive index (Fig. 4) as was proven for two well-known and popular explosives (RDX and NQ) and two materials seldom investigated in the terahertz domain (NTO and DNAN). Two unique absorption peaks at 1.54 and 2.75 THz can be observed for NTO, while DNAN is characterized only by one flat and relatively broad peak at 2.14 THz. Next, the transmission spectra of all IHEs were calculated with the use of characteristics of constituent materials and remained in good agreement with experimental measurements. This means that all peaks observed in the absorption characteristics have their corresponding slope on the refractive index curve, and moreover, all these features are directly connected with features of all constituents (Fig. 5). Of course, by comparing pure explosive materials and IHEs, some discrepancies can be distinguished and it is presumed that they result from local fluctuations in the concentration and thickness of materials.
The second method of terahertz examination of IHEs referred to the reflection spectra measurements. As before, we determined the characteristics of pure explosive materials (Fig. 6), but in comparison with the calculated spectra, their reflectance curves and phase characteristics were shifted towards higher frequencies. Summarizing, we can see that the reflectance characteristic corresponds to the refractive index variations and the reflectance phase to the absorption characteristic. Here, even a small misalignment of the sample surface or the reference mirror could cause some deviations from proper characteristics. The acceptable misalignment error is 10 μm.
Finally, the reflectance characteristics of IHE pellets were compared to the EMT-based calculations (Fig. 7). We think that the main role in the reflection measurements plays the surface-related scattering and not uniform surface distribution of the grains (Fig. 2a-c), which led to the proper reflectance characteristics being observed not at all measuring points. The surfaces of the IHE melt-castable discs were flat and uniform, and the terahertz radiation interacted with grains of RDX, NTO and NQ buried in the surface layer of the sample. In the case of RDX-based materials (PAX-41 and IMX-104), the solid solution of RDX in DNAN created the matrix, which enabled easier interaction with the particles. In the case of discs, we should also consider melt-related phenomena such as local concentration of the particles, clustering, breaking up of the particles and scattering at higher frequencies. The proposed effective medium theory offers a simple yet qualitatively correct description of the reflection from the complicated granular media. Better understanding of the related phenomena also requires consideration of the influence of the size of the terahertz beam. For 0.8 THz (RDX peak), the 1/e beam diameter is about 2 mm, while for 2.5 THz (NQ peak), it is below 0.5 mm. This factor can also explain the lower hit ratio of PAX-101 than for RDX-based IHEs.
The correlation coefficient between the measured and the theoretically calculated spectra has very high value for the transmission configuration. In all transmission cases, C is larger than 0.8 (Fig. 5) which is the sufficient confirmation of adequate choosing the EMT analysis method. For reflection measurements in most of the considered cases, the value of the correlation coefficient is high, but there are also some cases when it drops to small values (like for example in Figs. 6c and 7c), which results from the fact that the range 0.1-3 THz was used for correlation calculations. However, for both these cases, the similarity of the spectra is clearly seen, as well as the location of absorption peaks. For frequencies higher than 2 THz, often there are significant distortions of the spectrum which introduce a lower correlation value. Moreover, limiting the spectral range to 2 THz would mean that some distinctive spectral features would not be taken into account, which additionally depreciates the reliability of performed calculations. To compare such spectra, a more intelligent method as the principle component analysis can be used.
To sum up, we presented a comprehensive analysis of three melt-castable explosive materials and their constituents both in transmission and reflection. We proved that based on the refraction index of four constituent materials (DNAN, RDX, NTO and NQ) and EMT, we can calculate the reflection spectra of the IHEs with good agreement. In reflection, thanks to SEM images of the samples (pressed pellets and melt-castable discs), we explained the modifications of their reflectivity. In our opinion, grain-related scattering and clustering of particles must be considered. The presented analysis reveals that terahertz spectroscopy can be applied for reflection identification of IHEs. However, due to the inhomogeneity of the samples, measurement at a few (or even at a dozen) points must be considered to obtain reliable results. Moreover, since IHE materials are insensitive to compression, one can also measure their spectra by using an attenuated total reflection technique [27], in which sample preparation is easier than with pressed pellets.
In the future, the terahertz radiation can be used as a safe, fast and non-destructive method for quality control and examination of explosives. This radiation is absolutely non-ionizing and there is no risk of accidental initiation of the explosive being under tests. The results obtained by using the TDS method without doubt encourage further development of terahertz examination techniques and applying them to non-destructive analysis of explosives having a similar qualitative composition as the materials described in this paper. Moreover, this method allows for the detection of explosives containing relative compounds. Thus, the next step is to move from DNAN evaluation to analysis of trinitroanisole (TNAN), which can also be used as the melting ingredient.


extfieldsofstudy: Chemistry Medicine

text:
Luminescent Gold(III) Thiolates: Supramolecular Interactions Trigger and Control Switchable Photoemissions from Bimolecular Excited States

Abstract A new family of cyclometallated gold(III) thiolato complexes based on pyrazine‐centred pincer ligands has been prepared, (C^Npz^C)AuSR, where C^Npz^C=2,6‐bis(4‐ButC6H4)pyrazine dianion and R=Ph (1), C6H4 tBu‐4 (2), 2‐pyridyl (3), 1‐naphthyl (1‐Np, 4), 2‐Np (5), quinolinyl (Quin, 6), 4‐methylcoumarinyl (Coum, 7) and 1‐adamantyl (8). The complexes were isolated as yellow to red solids in high yields using mild synthetic conditions. The single‐crystal X‐ray structures revealed that the colour of the deep‐red solids is associated with the formation of a particular type of short (3.2–3.3 Å) intermolecular pyrazine⋅⋅⋅pyrazine π‐interactions. In some cases, yellow and red crystal polymorphs were formed; only the latter were emissive at room temperature. Combined NMR and UV/Vis techniques showed that the supramolecular π‐stacking interactions persist in solution and give rise to intense deep‐red photoluminescence. Monomeric molecules show vibronically structured green emissions at low temperature, assigned to ligand‐based 3IL(C^N^C) triplet emissions. By contrast, the unstructured red emissions correlate mainly with a 3LLCT(SR→{(C^Npz^C)2}) charge transfer transition from the thiolate ligand to the π⋅⋅⋅π dimerized pyrazine. Unusually, the π‐interactions can be influenced by sample treatment in solution, such that the emissions can switch reversibly from red to green. To our knowledge this is the first report of aggregation‐enhanced emission in gold(III) chemistry.

Synthesis and Characterization
When required, manipulations were performed using standard Schlenk techniques under dry nitrogen or in a MBraun glove box. Nitrogen was dried by passing through columns of supported P 2 O 5 , with moisture indicator, and activated 4 Å molecular sieves. Anhydrous solvents were freshly distilled from appropriate drying agents. (C^N pz^C )AuCl and (C^N py^C )AuCl were synthesised using literature methods. S1 1 H and 13 C{ 1 H} spectra were measured on a Bruker Avance DPX-300 spectrometer using CD 2 Cl 2 as solvent and signals were referenced to the residual protons of the deuterated solvent. IR spectra were recorded using a Perkin Elmer Spectrum 65 FT-IR spectrometer with a diamond ATR attachment. Elemental analyses were carried out at London Metropolitan University with a Perkin-Elmer 2400 CHNS/O microanalyzer.

Bu-4 (2)
Method 2: Under a nitrogen atmosphere, a flask was charged with (C^N pz^C )AuCl (56.0 mg, 0.097 mmol), K 2 CO 3 (26.9 mg, 0.191 mmol) and degassed acetone (5 mL). 4-tert-butylbenzenethiol (15 µl, 0.097 mmol) was added and the mixture was stirred at room temperature for 3 hours. The acetone was removed in vacuo and the resulting solid was extracted into dichloromethane (5 mL). The suspension was filtered through celite and washed with dichloromethane (5 mL). The solvent was removed. The solid residue was washed with cold hexane to give 2 as a yellow solid (56.7 mg, 0.081 mmol, 83%).

Synthesis of (C^N pz^C )AuSC 5 H 4 N-2 (3)
Method 1: (C^N pz^C )AuCl (50.0 mg, 0.087 mmol) and K 2 CO 3 (24.0 mg, 0.174 mmol) were dissolved in acetone (5 mL). Then, 1-mercaptopyridine (9.8 mg, 0.087 mmol) was added and the mixture was stirred at room temperature for 3 h. The acetone was removed in vacuo and the resulting solid was extracted into dichloromethane (5 mL). The suspension was filtered through celite and washed with dichloromethane (5 mL). The solvent was removed to give a yellow solid. This was washed with cold hexane to give 3 as a yellow powder (41.0 mg, 0.063 mmol, 72%). Single crystals suitable for X-ray diffraction were obtained from slow evaporation of a dichloromethane and isopropanol (4:1) solution.

Synthesis of (C^N pz^C )AuSAd (8)
Under an atmosphere of nitrogen, a flask was charged with (C^N pz^C )AuCl (55.0 mg, 0.097 mmol), K 2 CO 3 (26.9 mg, 0.191 mmol) and 1-adamantanethiol (16 mg, 0.097 mmol). Degassed acetone (5 mL) was added and the mixture was stirred at room temperature for 24 hours (the reaction was incomplete after 6 h, so an excess 1-adamantanethiol was added and the reaction was left overnight).

Aggregation of (C^N Pz^C )AuX complexes in solution
where I is the intensity of the observed spin echo, I 0 the intensity of the spin echo in the absence of gradient, D t the self-diffusion coefficient, Δ the delay between the midpoints of the gradients, δ the length of the gradient pulse, and γ the magnetogyric ratio. The shape of the gradients was rectangular, their length δ was 4-5 ms, and their strength G was varied during the experiments. All spectra were acquired for 64K points and a spectral width of 6200 Hz and processed with a line broadening of 1.0.
The semi-logarithmic plots of ln(I/I 0 ) versus G 2 were fitted by using a standard linear regression algorithm. Different values of G and number of transients were used for different samples.
The self-diffusion coefficient D t , which is directly proportional to the slope of the regression line where k is the Boltzmann constant, T is the temperature,  is the solution viscosity, c is the "size factor", which depends on the solute-solvent radius ratio and f is the "shape factor", which takes into account the deviation from sphericity of the diffusing molecule. a, b and d are the semi-axes of the ellipsoid that better approximates the shape of the diffusing molecule. S8 Structural parameter values (P= fc abd 3 = kT / phD t ) were derived from the measured D t values at different concentration in CD 2 Cl 2 and the data are collected in Table S3.1 and Figure S22.
Assuming a monomer/dimer equilibrium, the observed P values can be treated as the average between the parameters of both monomer (P M ) and dimer (P D ), according to the following relationship: where x M and x D are the molar fractions of monomer and dimer at a given concentration. P M and P D parameters are depending on the dimensions of monomer (a M , b M , d M ), dimer (a D , b D , d D ) and solvent (r solv , approximated at the van der Waals radius of CD 2 Cl 2 ). P M parameters can be derived from X-ray structures of monomers, by approximating their shape to ellipsoids. P D parameters have been calculated by using packing data derived by X ray diffraction (see below).  Figure S23. Evolution of the aromatic region of the 1 H NMR spectrum of (C^N Pz^C )AuSNp-2 as a function of the concentration in CD 2 Cl 2 .
The amount of dimers was quantified by using standard equations for monomer-dimer equilibrium: S9 where f m and f d are the molar fractions of monomers and dimers,  M and  D are the limit chemical shift values of monomer and dimer, respectively, K is the association constant. After fitting chemical shift versus concentration trends, the molar fractions of dimers were obtained as: Table S3.2.

Photophysical Properties
UV−visible absorption spectra were recorded using a Perkin-Elmer Lambda 35 UV/vis spectrometer.
Excitation and emission spectra were measured using a (

Cryogenic time resolved photoluminescence (ICCD and TCSPC)
ICCD : An electrically gated intensified CCD ( ICCD ) camera (Andor iStar DH740 CCI-010) with a calibrated grating spectrometer (Andor SR303i) were used for the time resolved PL spectra. 400nm laser light pulsed at 1 KHz was used to pump the sample. This was produced from the second harmonic generation of a 800nm femtosecond laser using a BBO crystal. The fundamental output had a 80 fs pulse width and 1.5 eV pule energy from a Ti:Sapphire laser system (Spectra Physics Solstice).
Time resolution of the emission was achieved by changing the delay of the gate relative to the 1 KHz pulsed pump laser. The minimum gating width was ~2.5 ns. A 495 nm long pass filter was used to remove any scattered light from the pump laser. This setup was used for the green emitting film.
Time correlated single photon counting : A pulsed 407 nm laser was used to photo-excite the samples, with pulse width <200 ps and a rep rate of 2.5MHz. The detector was a Si based single photon avalanche photodiode. The instrument response function was ~200 ps also. A 495 nm long pass filter was used as in the ICCD. This setup was used for the red emitting film.

Theoretical Calculations
All structures were optimised using density functional theory, within the approximation of the PBE0 exchange and correlation functional S11 as implemented within the ORCA quantum chemistry program. S12 A def2-SVP basis set S13 was used throughout. To account for the weak π−π interactions, the calculations were supplemented with Grimme's D3 dispersion correction with the Becke−Johnson damping scheme. S14 The absorption spectrum of complex 5 was calculated using pertubrative spin orbit coupling TDDFT (pSOC-TDDFT) S15 calculated using the ZORA Hamiltonian S16 as implemented within ADF. S17 The PBE0 exchange and correlation functional was used, and a DZP basis set for all atoms, except gold, for which a TZP basis was used. All emission energies were calculated using both uDFT and pSOC-TDDFT within ADF.


text:
A new anti-inflammatory β-carboline alkaloid from the hairy-root cultures of Eurycoma longifolia

Abstract One new β-carboline alkaloid 7-methoxy-(9H-β-carbolin-1-il)-(E)-1-propenoic acid (1) together with 9-methoxycanthin-6-one (2) and 9-hydroxycanthin-6-one (3) were isolated from the hairy-root cultures of Eurycoma longifolia. The effects of these compounds on lipopolysaccharide (LPS)-induced nitric oxide (NO) production in RAW264.7 cells were investigated. Compound 1 strongly inhibited the production of NO while 2 and 3 having weak or inactive effect. Consistently, compound 1 decreased the expression of cyclooxygenase-2 and inducible nitric oxide synthase. Graphical abstract

Introduction
Eurycoma longifolia Jack (Simaroubaceae) is an important medicinal plant in South-East Asian countries due to its broad pharmaceutical properties such as anti-inflammation, anticancer, antidiabetes and antimalaria (Bhat & Karim 2010). The roots of this plant have been widely used as a single ingredient or as a part of a mixture with other herbs for the treatment of sexual insufficiency, malaria, aches, dysentery and fever. In Vietnam, E. longifolia is named 'cay ba benh' , meaning a tree that cures hundreds of diseases (Loi 2004). Previous studies have reported that the plant contains quassinoids, canthin-6-one and β-carboline alkaloids, squalene derivatives, tirucallane-type triterpenes and biphenylneolignans (Kuo et al. 2003;Park et al. 2014;Tran et al. 2014). However, E. longifolia roots are harvested after several years of cultivation, i.e. time-consuming and seasonal depending, and the uprooting harvest for medicinal preparation has caused an adverse effect on the natural population of E. longifolia. Thus, an alternative method of E. longifolia production is required to meet the ever-growing demand of this medicine. To solve this problem, the use of plant tissue culture to produce high-value secondary compounds for pharmaceutical and industrial purposes is ever increasing. In comparison with different plant tissues, hairy-root cultures have been reported many advantages such as continuous growth periods, rapid growth rates, easy culturing in hormone-free media and accessibility to genetic manipulation and thus is a useful approach for the production of active compound in confined and secure environments. Hairy roots can be obtained by infecting the wounded higher plants with Agrobacterium rhizogenes, a gram negative soil bacterium that can transfer a DNA segment (T-DNA) from root-inducing (Ri) plasmid into the genome of the infected plant. This T-DNA carries a set of genes that encode enzymes which control auxin and cytokinin biosynthesis, which induces the formation of hairy roots (Guillon et al. 2006). In this study, cultured hairy roots of E. longifolia were produced using the A. rhizogenes-mediated transformation system. A chromatographic fractionation of 90 g dried mass of hairy roots led to the isolation of a new β-carboline alkaloid (1) together with 9-methoxycanthin-6-one (2) and 9-hydroxycanthin-6-one (3) (Kardono et al. 1991).

Results and discussion
Compound 1 was obtained as amorphous yellow powder. Its HR-ESI-MS spectrum revealed the peak at m/z 269.0929 [M+H] + corresponding to the molecular formula C 15 H 13 N 2 O 3 . The 1 H-NMR spectrum of 1 showed the signals characteristic of an ABX aromatic system at δ H 6.95 ( . The 13 C-NMR and DEPT spectra of 1 indicated the presence of a methoxy, seven methines and seven quaternary carbons. A carboxylic (δ C 166.8) and two methine (δ C 128.9 and 129.1) signals were attributed for a propenoic acid moiety. The rest of signals corresponded to the β-carboline skeleton previously found in Eurycoma sp. (Kuo et al. 2003;Park et al. 2014;Tran et al. 2014). These data were very similar to those of 6-methoxy-(9H-β-carbolin-1-il)-(Z)-2-propenoic acid as previously reported (Nunomura et al. 2012). However, the large coupling constant (J = 13.0 Hz) of the olefinic protons H-1′ and 2′ indicated the trans double bond in 1. Extensive analysis of HMBC correlations ( Figure 1) confirmed that the methoxy group attached to C-7. Thus, 1 was elucidated to be 7-methoxy-(9H-β-carbolin-1-il)-(E)-1-propenoic acid.
All of the isolated compounds were tested for their inhibition of NO production in LPSstimulated RAW264.7 cells. As a result, the β-carboline 1 strongly prevented the NO production in a dose-dependent manner while two canthin-6-ones 2 and 3 were weak or inactive up to 30 μM ( Figure 2). The cell viability assay indicated that the inhibition of NO production was not resulted from cytotoxic effects of the tested compounds (data not shown). It is known that cyclooxygenase-2 (COX-2) and inducible nitric oxide synthase (iNOS) are two major inflammatory mediators, and the induced expression of iNOS is primarily responsible for the NO production. Thus, the effects of 1 on LPS-induced expressions of COX-2 and iNOS in RAW264.7 cells were also investigated using Western blot analysis. As the result, 1 dose-dependently suppressed the LPS-induced expressions of iNOS and COX-2 proteins. It has been reported that β-carboline alkaloids exert extensive biochemical activities and multiple pharmacological effects such as anti-inflammation, antitumor, antibacteria and antithrombosis (Lee et al. 2000;Cao et al. 2007;Zhang & Sun 2015). Consistently, our results showed that the β-carboline 7-methoxy-(9H-β-carbolin-1-il)-(E)-1-propenoic acid (1) strongly down-regulated the inflammatory iNOS and COX-2 expression.
In conclusion, this work identifies a new anti-inflammatory β-carboline, 7-methoxy-(9Hβ-carbolin-1-il)-(E)-1-propenoic acid, from the cultured hairy roots of E. longifolia. This compound showed strong inhibitory effect on the NO production and decreased COX-2 and iNOS protein expression in LPS-induced RAW264.7 cells. Further studies are in progress to determine the mode of action of this anti-inflammatory agent.

General procedures
TLC was performed using precoated Kiesel gel 60 F254 and visualised by UV light 254 nm and 10% H 2 SO 4 reagent with heat. Column chromatography was performed using silica gel 60 (Merck, 70-230 mesh). NMR experiments were carried out on a Bruker AM500 FT-NMR spectrometer using tetramethylsilane as internal standard. The HR-ESI-MS spectra were recorded on an API Q-STAR PULSAR I of Applied Biosystem.

Plant materials
The seeds of E. longifolia Jack were collected in July 2011 in Bai Tu Long National Park, Quang Ninh, Vietnam, and identified by Dr Bui Van Thanh, Institute of Ecology and Biological Resources, Vietnam Academy of Science and Technology. The voucher specimens (BTL11) were deposited at the Institute of Biotechnology.

Preparation of cultured hairy roots
Seeds of E. longifolia Jack were obtained from a mature tree in Bai Tu Long National Park, Quang Ninh province, Vietnam. Seeds were washed thoroughly with detergent and rinsed under running tap water for 5 min. After removing the seed coat, the seeds were sterilised with 70% ethanol for 1 min and 0.1% mercuric chloride solution for 3 min, followed by rinsing with autoclaved water. Sterilised seeds were germinated on MS medium supplemented with 3% sucrose and GA3 (0.3 mg/L). Cotyledons and hypocotyls from 3-week-old in vitro cultures were used for the A. rhizogenes transformation.
A colony of A. rhizogenes strain ATTC 15834 (provided by Institute of Pharmacy and Molecular Biotechnology, Heidelberg University, Germany) was cultured on 50 mL YMB medium at 28 °C, 200 rpm on an orbital shaker for 14-16 h. Then, the bacterial cells were collected by centrifugation at 5000 rpm for 10 min and resuspended in 25 mL ½ Woody Plant (WP) medium. Cotyledons and hypocotyls were cut in 0.3-0.4 cm in ½ WP liquid medium. Explants dried on sterile filter papers were transferred onto a new plate containing 25 mL of agrobacterial suspension and incubated for 30 min with gentle shaking. The bacteria were poured off, and the explants were dried on sterile filter paper and placed on WP solid medium for 2 days under low light intensity. After three days inoculated with bacteria, the explants were placed onto new WP plates containing 500 mg/L cefotaxime (Brithol Michcoma, Holland). The explants were subcultured on new WP antibiotic plates every two weeks until hairy root tips appeared. Hairy root tips (4-5 cm in length) were then separately transferred to 50-mL flasks containing 10 mL of WP liquid medium and placed on a shaker at 60 rpm, 25 ± 2 °C for further analysis. A hairy root exerting high growth rate was selected to produce hairy root biomass. Four grams of hairy roots biomass were cultured in 150 mL liquid WP medium (semi-submerged culture) at 25 °C and orbitally shaked 60-70 rpm in the dark. After 3 weeks of cultures, 200 mL of WP fresh medium was supplied into the cultures. After 6 cultured weeks, the biomass obtained was dried at 37 °C and average production of hairy root (5 g dry weight, DW). Hairy root was continuously cultured for high amount of mass, and the dry mass was used for next experiments.

Extraction and isolation
The cultured hairy roots of E. longifolia (90 g dried mass) were extracted with methanol in a sonic bath (500 mL × 3 times). The combined extracts were concentrated to give 13.5 g crude extract which was then partitioned between 10% HCl solution and ethyl acetate (each 500 mL). After the removal of organic layer, the acidic solution was adjusted to pH 10 by 1 N NaOH and then extracted with ethyl acetate (300 mL × 3 times) to obtain an alkaloid residue (651 mg). The alkaloid extract was chromatographed on a silica gel column eluted with dichloromethane-methanol (25:1 v/v) to afford 12 fractions F1-12. Compound 3 (15.0 mg) was recrystallised in dichloromethane from F4. Compound 2 (96.0 mg) was purified from F2 by a C18 column using methanol-water 3:1 (v/v). Fraction F12 was chromatographed on a silica gel column eluted with dichloromethane-methanol (20:1 v/v) to obtain 1 (29.6 mg).

Assay for inhibition of NO production and iNOS, COX-2 expression
Determination of NO production and the immunoblot analysis for iNOS and COX-2 were performed as previously described (Trang et al. 2014). Celastrol was used as a positive control for the NO assay (IC 50 , 1.0 μM).


extfieldsofstudy:
Materials Science
Medicine

text:
Cooperative Catalytic Effect of ZrO2 and a-Fe2O3 Nanoparticles on BiVO4 Photoanodes for Enhanced Photoelectrochemical Water Splitting

at present. To overcome this problem, different strategies have been implemented, including ion doping, g, 4a, c, 5] nanostructuring, f, 4a, 6] surface modification with passivation layers or electrocatalysts, d, 5d, 7] and combinatorial synthesis. For example, a-Fe2O3 photoanodes have been successfully doped with heteroatoms like Sn, Ti, Zr, c] Si, and Nb. In the case of monoclinic BiVO4 (m-BiVO4), &&definition added.&& there are two different sites susceptible to metal doping, that is, Bi + and V + . Whereas the replacement of V sites by six valence metal ions, such as Mo and W, 5e, f] has been largely explored, there are very few studies on doping Bi + sites, and are generally limited to theoretical calculations. Hf + , Zr , and Sn + , with ion radius relatively similar to Bi , are potential candidates to occupy Bi positions in the m-BiVO4 crystal structure. Luo et al. synthesized BiVO4 photoanodes doped by Sn 4 in Bi 3 positions by means of a metal–organic decompo-

Photoelectrochemical water splitting with metal oxide semiconductors offers a cost-competitive alternative for the generation of solar fuels. Most of the materials studied so far suffer from poor charge-transfer kinetics at the semiconductor/liquid interface, making compulsory the use of catalytic layers to overcome the large overpotentials required for the water oxidation reaction. Herein, we report a very soft electrolytic synthesis deposition method, which allows remarkably enhanced water oxidation kinetics of BiVO 4 photoanodes by the sequential addition of Zr and Fe precursors. Upon a heat treatment cycle, these precursors are converted into monoclinic ZrO 2 and a-Fe 2 O 3 nanoparticles, which mainly act as catalysts, leading to a five-fold increase of the water oxidation photocurrent of BiVO 4 . This method provides a versatile platform that is easy to apply to different semiconductor materials, fully reproducible, and facile to scale-up on large area conductive substrates with attractive implications for technological deployment.
&&Minor change to title for purely formatting reasons. Change ok?&& Photoelectrochemical (PEC) water splitting driven by sunlight constitutes an environmentally friendly technology to store energy in the covalent bonds of H 2 . [1] To efficiently convert solar photons into such a chemical fuel, semiconducting materials capable of absorbing a large fraction of the solar spectrum, with an adequate band edge position, and low overpotentials to carry out the hydrogen and/or oxygen evolution reactions are required. In particular, intensive research efforts have been accomplished to develop competitive n-type semiconducting materials working as water splitting photoanodes, as the oxygen evolution reaction (OER) is both kinetically and thermodynamically more demanding, [2] and the material suffers from severe oxidizing conditions. To achieve cost-competitive, efficient, and durable devices, the use of low-cost, earth-abundant, stable materials synthesized by easily up-scalable methods is essential. In this context, metal oxides, such as monoclinic BiVO 4 (m-BiVO 4 ) [3] and a-Fe 2 O 3 , [4] have been deeply investigated as candidate photoanode materials. Although both of them satisfy many of the requirements, their poor electronic properties yield low solar-tofuel efficiencies, preventing their large-scale use in PEC systems at present. To overcome this problem, different strategies have been implemented, including ion doping, [3d, g, 4a, c, 5] nanostructuring, [3a, f, 4a, 6] surface modification with passivation layers or electrocatalysts, [3b, d, 5d, 7] and combinatorial synthesis. [8] For example, a-Fe 2 O 3 photoanodes have been successfully doped with heteroatoms like Sn, [5a] Ti, [4c] Zr, [5b, c] Si, [4a] and Nb. [5d] In the case of monoclinic BiVO 4 (m-BiVO 4 ), &&definition added.&& there are two different sites susceptible to metal doping, that is, Bi 3 + and V 5 + . Whereas the replacement of V 5 + sites by six valence metal ions, such as Mo and W, [3d, 5e, f] has been largely explored, there are very few studies on doping Bi 3 + sites, and are generally limited to theoretical calculations. [9] Hf 4 + , Zr 4 + , and Sn 4 + , with ion radius relatively similar to Bi + 3 , are potential candidates to occupy Bi positions in the m-BiVO 4 crystal structure. Luo et al. [10] synthesized BiVO 4 photoanodes doped by Sn + 4 in Bi + 3 positions by means of a metal-organic decomposition (MOD) method. A practically negligible improvement in photocurrent was reported for the doped sample, which was attributed to the high formation energy and the relatively low solubility of Sn 4 + impurity ions into the crystal lattice.
The functional properties of BiVO 4 and Fe 2 O 3 have also been separately improved by means of heterostructuring strategies. For instance, the PEC performance of BiVO 4 has been largely increased when deposited on nanostructured WO 3 layers, owing to the synergistic interaction between the BiVO 4 (providing good light harvesting properties) and WO 3 (enhanced charge transport). [3f, 11] Other different host-guest semiconductor combinations have been successfully implemented and tested with remarkable results for water splitting, like Si/TiO 2 / BiVO 4 , [12] SnO 2 /BiVO 4 , [13] WO 3 /Fe 2 O 3 , [14] and Si/Fe 2 O 3 . [15] Nevertheless, to the best of our knowledge, the heterostructured Fe 2 O 3 /BiVO 4 system, which would entirely match the stability and low-cost requirements, has not been tested for solar water oxidation. There is only one recent work in literature reporting the decoration of BiVO 4 with Fe 2 O 3 nanoparticles, where it is suggested that Fe 2 O 3 acts as an efficient co-catalyst for the degradation of organic pollutants. [16] Although different Fe-based catalysts have been successfully employed to deco- rate semiconductor photoanodes, such as FeOOH [7a] and Ni-FeO x , [3a, 17] the use of a-Fe 2 O 3 nanoparticles as OER catalyst has not been explored. Herein, we report a very soft electrolytic synthesis deposition method that is fully reproducible and facile to scale-up on large area conductive substrates, which allows enhancing remarkably the water oxidation kinetics of BiVO 4 photoanodes by the sequential addition of Zr and Fe precursors.
BiVO 4 photoanodes were synthesized by means of the twostep method developed by Choi et al., [3c] which consisted of Bi electrodeposition on a fluorine-doped tin oxide (FTO) &&definition added.&& substrate followed by drop casting &&ok? && the V precursor and annealing at 500 8C for 2 h. Different amounts of Zr 2 Cl 2 O were incorporated into the Bi 3 + plating bath as the source of Zr. Subsequently, electrodeposition of Fe on the as-prepared BiVO 4 layers followed by annealing at 450 8C was carried out to obtain Fe 2 O 3 nanoparticles. Different electrodeposition charges were tested, ranging from 1 to 32 mC cm À2 . After systematic optimization of the deposition process, a five-fold increase of the water oxidation photocurrent at 1.23 V versus RHE was obtained compared to the reference BiVO 4 system. Figure 1 illustrates the photoelectrochemical behavior of the reference Fe 2 O 3 (A) and BiVO 4 (B) photoanodes together with the best modified BiVO 4 electrodes: BiVO 4 -Zr (C), BiVO 4 -Fe (D), and BiVO 4 -Zr-Fe (E). The photocurrents at 1.23 V versus RHE for the different Zr and Fe additions are showed in Figure 1 b and are used &&ok?&& to determine the optimum synthetic conditions reported in Figure 1 a. Typical volcano plots for both Zr and Fe additions on top of the BiVO 4 surface were obtained. The best performance was obtained for 2.5 mole % Zr. Note that this concentration is referred to the Zr/Bi molar ratio added to the electrodeposition bath. On the other hand, 2 mC cm À2 of total charge deposition of Fe yielded the optimum results. Although the photocurrents reported in Figure 1 are lower compared to previous studies on BiVO 4 photoanodes, [3a, b, 7a] the efficiency is remarkably enhanced by the Zr and Fe additions compared with the reference material. Additionally, the simplicity of the method de-scribed here to engineer surface modifications can be mimicked in more efficient systems. Figure 1 c shows the j-V curves when a hole scavenger (1 m Na 2 SO 3 ) was added to the solution. In general, all BiVO 4 -based samples exhibit practically the same behavior for sulfite oxidation when surface recombination is negligible, with photocurrents slightly higher than reported for oxygen evolution. This strongly suggests that the improvement reported for the optimized BiVO 4 -Zr-Fe photoanode is mainly connected to surface catalysis.
To more precisely determine the role of both additives (Zr and Fe) on the PEC behavior of the photoelectrodes, detailed structural and optical characterization was performed. First, XRD measurements of the metallic deposits were carried out (see Figure SI1 in the Supporting Information). The results show that as the Zr content increases, the intensity of the peak assigned for Bi at 2q = 27.58 diminishes, which indicates that Bi-Zr codeposition was successfully accomplished. In addition, SEM images of as-deposited metals ( Figure SI2) show a progressive morphological evolution with increased Zr additions, from dendrite trunks and branches (pure Bi deposits) to more compact films. Figure 2 a and b show the SEM micrographs of pristine BiVO 4 and BiVO 4 -Zr, with the optimum Zr addition (2.5 mol %), illustrating that after the vanadate formation process, the microstructure of the samples considerably changes to smoother-edged dendrites, although a high surface area is maintained. Furthermore, it is remarkable the decrease of the nanoparticle size from approximately 500 nm for pristine BiVO 4 to 200 nm for the optimized BiVO 4 -Zr electrode. This is consistent with the role of Zr as grain growth inhibitor as reported for different metallurgical synthetic processes. [18] Figure 2 c shows the XRD diffractograms of BiVO 4 -Zr samples for different Zr additions, perfectly matching the monoclinic scheelite BiVO 4 structure (CAS Number 00-014-0688). The intensity of the (1 2 1) &&ok?&& peak increases with the Zr concentration up to the optimum 2.5 mol %, and is practically unchanged with further Zr additions. Moreover, the peak position shifts toward higher 2q values for the most Zr-rich samples. This clearly involves a degradation of the film integrity as  Communications   1  1  2  2  3  3  4  4  5  5  6  6  7  7  8  8  9  9  10  10  11  11  12  12  13  13  14  14  15  15  16  16  17  17  18  18  19  19  20  20  21  21  22  22  23  23  24  24  25  25  26  26  27  27  28  28  29  29  30  30  31  31  32  32  33  33  34  34  35  35  36  36  37  37  38  38  39  39  40  40  41  41  42  42  43  43  44  44  45  45  46  46  47  47  48  48  49  49  50  50  51  51  52  52  53  53  54  54  55  55  56 56 57 57 a result of the strain induced by the smaller ionic radius of Zr (0.79 ) replacing Bi 3 + positions (1.11 ) into the periodic crystal lattice, [5g] which can be related to the progressive decrease of the photocurrent showed in Figure 1 b for the higher Zr additions. Figure 3 shows HR-TEM images of the optimized BiVO 4 -Zr and BiVO 4 -Zr-Fe photoanodes. The BiVO 4 -Zr electrodes are characterized by the presence of high crystalline nanoparticles (5-10 nm), with interplanar distances of 2.8 (Figure 3 a) and 3.1 (Figure 3 b) corresponding to the (111) and (1 11) reflections of monoclinic ZrO 2 (CAS Number 1309-37-1), respectively. Local energy dispersive X-ray spectroscopy (EDS) &&definition added.&& analyses confirmed the presence of Zr at these locations ( Figure SI3), clearly indicating that although Zr can substitute Bi in the monoclinic BiVO 4 lattice as showed by XRD, an important fraction of Zr is present at the surface of BiVO 4 in the form of monoclinic ZrO 2 nanoparticles. The formation of ZrO 2 nanoparticles upon Zr addition has been reported even for low doping densities, below 1 at %. [19] On the other hand, the optimized BiVO 4 -Zr-Fe samples, additionally showed the presence of crystalline nanoparticles, with interplanar distances of 2.7 (Figure 3 c), corresponding to the (1 0 4) reflections of hexagonal scalenohedral a-Fe 2 O 3 , (CAS Number 1309-37-1) that appeared homogeneously dispersed on top of the BiVO 4 -Zr photoanode.
The optical properties of BiVO 4 were not significantly altered by Zr additions, as showed in the Supporting Information (Figure SI4 a), in good agreement with previous reports. [20] According to its bandgap of 2.4 eV, BiVO 4 is able to absorb light up to approximately 510 nm. On the other hand, as the amount of electrodeposited Fe increases, a broad shoulder appears on the 500-650 nm region, matching the main absorption band of a-Fe 2 O 3 (Figure SI4 Communications   1  1  2  2  3  3  4  4  5  5  6  6  7  7  8  8  9  9  10  10  11  11  12  12  13  13  14  14  15  15  16  16  17  17  18  18  19  19  20  20  21  21  22  22  23  23  24  24  25  25  26  26  27  27  28  28  29  29  30  30  31  31  32  32  33  33  34  34  35  35  36  36  37  37  38  38  39  39  40  40  41  41  42  42  43  43  44  44  45  45  46  46  47  47  48  48  49  49  50  50  51  51  52  52  53  53  54  54  55  55  56  56  57  57 agreement with those previously reported. A progressive decrease of the bandgap energy with the addition of Fe was detected, from an electrodeposition charge of 12 mC cm À2 , obtaining E g = 2.32 eV for the highest Fe loading tested (i.e., 32 mC cm À2 ), in good agreement with the optical properties of a-Fe 2 O 3 . From these optical absorption measurements, the maximum achievable photocurrent was calculated considering that all the photogenerated holes participate in the water oxidation reaction (j abs ), see Table SI1. The enhanced PEC activity induced by the Zr and Fe additions was further corroborated by the incident photon-to-current efficiency (IPCE) spectra obtained at 1.23 V versus RHE without and with the addition of the hole scavenger (Figure 4 a  and b, respectively). There is an excellent correspondence between the onset wavelength of the IPCE and the absorbance measurements. Moreover, the integrated photocurrents extracted from the IPCE spectra nicely match those obtained from the j-V curves, as summarized in Table SI1. From the IPCE measurements in the solution with the hole scavenger, it is clear that the enhanced PEC behavior for the optimized heterostructure can be ascribed to the "cooperative" (rather than synergistic) catalytic effect of both ZrO 2 and a-Fe 2 O 3 nanoparticles deposited on the top of the m-BiVO 4 films. Evidence of the excellent electrocatalytic behavior of a-Fe 2 O 3 has been previously reported for Si photoanodes. [21] Conversely, nanostruc-tured monoclinic ZrO 2 nanoparticles have demonstrated catalytic activity for different chemical reactions, [22] and has been also employed to passivate TiO 2 surface traps in dye-sensitized solar cells. [23] However, to the best of the authors knowledge it is the first time that they show enhanced water oxidation kinetics.
To quantitatively assess the beneficial effect of both Zr and Fe additions on BiVO 4 films, the charge separation (h cs ) and charge injection (h cat ) yields were calculated for the synthetized materials through Equations (1) and (2), by comparing the photocurrent for water splitting and to a hole scavenger: [24] The comparative behavior of both calculated yields is showed in Figure 5. When Fe is distributed in the form of scattered Fe 2 O 3 nanoparticles on the BiVO 4 surface, this material acts as an excellent catalyst for oxygen evolution reaction, doubling the catalytic efficiency of the semiconductor substrate. It is also remarkable the increase in the catalytic efficiency of BiVO 4 -Zr samples, with respect to the reference BiVO 4 . In contrast, the effect of Zr addition on the charge separation effi-   Communications   1  1  2  2  3  3  4  4  5  5  6  6  7  7  8  8  9  9  10  10  11  11  12  12  13  13  14  14  15  15  16  16  17  17  18  18  19  19  20  20  21  21  22  22  23  23  24  24  25  25  26  26  27  27  28  28  29  29  30  30  31  31  32  32  33  33  34  34  35  35  36  36  37  37  38  38  39  39  40  40  41  41  42  42  43  43  44  44  45  45  46  46  47  47  48  48  49  49  50  50  51  51  52  52  53  53  54  54  55  55  56  56  57  57 ciency is negligible, in good agreement with previous reports suggesting that Zr additions on BiVO 4 do not affect the electronic properties of the material. [3e, 20] The validity of this analysis is based on the assumption of complete absence of surface recombination when the hole scavenger is employed. A marginal difference is observed in the measurements with the hole scavenger, which can be attributed to the role of ZrO 2 and a-Fe 2 O 3 as surface passivation layers on top of BiVO 4 . [7b] In any case, further analysis is needed to clarify this issue. Figure SI6 shows the Mott-Schotky plots obtained in the dark on BiVO 4 and optimized BiVO 4 -Zr and BiVO 4 -Zr-Fe samples. The lower slope reported for the BiVO 4 -Zr samples is connected to an increase in the surface area as shown in Figure 2, rather than to an increase in donor density. The material roughness is further increased for the BiVO 4 -Zr-Fe sample, which can explain the lowest slope in the Mott-Schottky plot of Figure SI6, as a consequence of the Fe 2 O 3 nanoparticle decoration and/or the ZrO 2 morphological rearrangement on top of BiVO 4 photoanode under the second annealing. [25] The calculated flat band potentials (V fb ) and donor densities (N D ) are collected in Table SI2. No significant changes are obtained for the electronic properties of BiVO 4 upon Zr and Fe additions, as reflected in the charge separation yield in Figure 5 a.
To benchmark the cooperative catalytic Zr-Fe-based layer developed in the present study, we have compared the reported behavior with that obtained with a FeOOH catalyst. For this purpose, we have deposited a FeOOH layer on top of our BiVO 4 photoanode, according to the optimum synthesis conditions described by Seabold and Choi in Ref. [7a]. Figure SI7 shows both the photoelectrochemical and optical properties of both systems. Our BiVO 4 -Zr-Fe photoanode clearly outperforms BiVO 4 -FeOOH in terms of delivered photocurrent. Additionally, the absorbance of our catalytic layer is significantly lower, which is particularly beneficial for application in tandem architectures.
In summary, we have reported a facile method to enhance the PEC water oxidation behavior of BiVO 4 photoanodes by controlled additions of Zr and Fe precursors on the electrodeposition bath. A remarkable five-fold increase of the photocurrent is reported for the optimized BiVO 4 -Zr-Fe photoanodes, which can be explained by the cooperative catalytic role of monoclinic ZrO 2 and a-Fe 2 O 3 nanoparticles distributed on the surface of BiVO 4 . Although we show by XRD analysis that Zr can also extrinsically dope BiVO 4 replacing Bi atoms, the effect on the intrinsic electronic properties of BiVO 4 is not beneficial for the optimum Zr addition. The findings of this study can also be easily extrapolated to enhance the PEC performance of different photoanode materials with low-cost catalytic materials like Fe 2 O 3 and ZrO 2 , obtained by a simple, fully reproducible, and up-scalable synthesis.