QM/MM methods in studies of coinage metals: copper, silver, and gold interacting with biological and organic molecules

ABSTRACT

A QM/MM method is an atomistic simulation algorithm that allows researchers to describe different regions of a system with different physical laws. Here, we review this hybrid method’s applications to the study of copper, silver, and gold atoms and clusters interacting with biological and organic molecules. In particular, we highlight efforts to characterize the relaxation process in a copper(I) phenanthroline complex; details of Cu’s secretory path; the atomic structure of Ag-homopolymers of cytosine and guanine; DNA-stabilized silver clusters; effects related to temperature and solvent on thiolate-protected gold clusters’ optical properties; and the effect of a medium-like noble gas on a cluster’s optical spectrum. The results of these efforts demonstrate how QM/MM methods are applied successfully to a wide range of processes that include the study of excited state evolution, charge transport, light absorption, and emission, and determining an atomic structure in the absence of crystal-determined structure. We expect QM/MM methods will continue supporting the exploration of novel hybrid organo-metallic materials and their safe use in the environment, while also providing guidance on mechanisms to deal with diseases associated with a failure in cells’ proper behavior.

Graphical Abstract

KEYWORDS:

• QM/MM
• metal-organics
• Copper-mediated
• DNA-stabilized clusters
• thiolate-protected gold clusters

1. Introduction

Metals exist in many fundamental biological processes and form hybrid complexes with organic matter. Copper, for example, is essential to copper-mediated enzymes that regulate oxygen transport and communication between neurons [1], and iron is present in hemoglobin and oxygen transport processes [2]. To date, nanoscience has offered some exploration of hybrid metalo-organic materials, such as DNA-stabilized silver clusters [3–5] and thiolate-protected gold clusters [6,7] with potential applications as biosensors, in bioimaging and medicine [8,9]. Yet within these systems, many properties and possibilities remain unexplored. Other knowledge we might gain about these metals in hybrid form include questions on their stability in solution, which are important questions in potentiating their use in technological applications.

To better understand and push beyond the body of research amassed thus far, here we focus on hybrid materials made of Cu, Ag, and Au metals and biological or organic matter. The metal atoms are part of the chemical periodic table’s group 11 and also belong to the noble and coinage metals category. As part of the same column, the metal atoms share a common pattern in their electronic structure. The outermost shell is confirmed by a single s-type electron, and the shell energetically below is a d-type closed shell. When such atoms aggregate in the nano regime–i.e. in atomic combinations where the scale length is nanometric–the electron–electron interaction gives rise to a mean-field type of model termed the shell model [10]. In this configuration, the 1s electron forms a delocalized band, but because of the model’s quantum nature, it has defined shells that are filled by such electrons. In the shell model, therefore, the number of atoms or electrons in the s-band determine the filling of the shells. Nanosystems with closed shells have increased stability. This increased stability at a given number of atoms is termed ‘magic sizes.’ A consequence of this type of interaction is that it needs a quantum description.

Typical biomolecules tend to have large-scale structures with a wide space of conformations that follow dynamical processes via electrostatic interactions, and metals participate with their reactive electronic structure where quantum effects dominate. Depending on the ligand’s chemical nature, the interaction between the metal and organic matter can be covalent (for example, the Au-S bond in thiolate-protected clusters with shared electrons forming a chemical bond) or ionic (where a full electron transfers typically from the 1s outermost shell to the ligand, with a typical example being the Ag-N bond in Ag-DNA helices). This type of interaction modifies the electronic properties, such as optical absorption and the circular dichroism spectrum. When metals are in solution, the polarization interaction might lead to important effects. In general, the hybrid nature of these systems makes their complex chemical processes to be unwieldy to determine, but if we use existing hybrid methods for computational simulations, this becomes an approachable challenge.

Computational simulations are one of the most powerful tools for studying chemical processes. They provide atomistic detail of complex mechanisms and improve insights into what can be gained from experimental results. However, computational results have limitations in practice. In principle, quantum mechanics is the theory that approximates molecules and nanomaterials most accurately. A fast and accurate enough quantum method such as density functional theory (DFT) [11] can only reach on the order of hundreds of atoms in a timescale of picoseconds simulations. This clearly prevents applying this methodology to systems such as large biological molecules.

On the other hand, a cheaper method to simulate systems of a larger scale regarding the number of atoms and time is the classical approximation, which uses approximated force fields for modeling the interatomic forces that lead the system’s changes. The principal problem with classical approximation is that it improves the simulated scales but at the cost of precision in results, mainly in systems with metals where electronic s and d electrons need a quantum description.

A hybrid method called quantum mechanics/molecular mechanics (QM/MM) is thus the method of choice to perform simulations where some part of the system will be described at the quantum level and another in the classical level. A wide range of codes implement this method, including SIESTA [12], CHARMM [13], GPAW [14], CPMD [15], CP2K [16] and AMBER [17]. Current developments to improve QM/MM methods seek the use of semi-empirical quantum methods or energies and forces from neural networks trained on quantum simulations [18].

QM/MM has been used to investigate several metal-organic compounds, and although we focus here in coinage metals, we want to mention that there exist applications of this method to other metals, among them, Ru [19–23], Pd [24–27] and Pt [28–31]. An example of an implementation of QM/MM to research Ru metal-organic compounds is the study of the Ru-tris 2,2’-bipyridine in excited states that undergo to metal to ligand charge transfer (MLCT) state [22,23]. Another example is the investigation of the kinetic of [(-benzene)Ru(II)(en)(H O)] (en = ethylenediamine) binding mechanism to the N7 position of guanine and a subsequent cross-Link formation between adjacent guanines producing DNA deformation [20]. For the case of Pd, QM/MM has been used to address the solvation of nanoparticles in ionic liquids [26,27].

Here, we will review studies that characterize the physicochemical atomic and electronic properties of metal-organic Cu, Ag, and Au, including both atomic and cluster metal forms. We highlight efforts to characterize the relaxation process in a copper(I) phenanthroline complex; Ag-homopolymers of cytosine and guanine; DNA-stabilized fluorescent silver clusters; thiolate-protected gold clusters’ and noble gas-stabilized silver clusters. To further the exploration of these metals’ hybrid properties and potential, here we provide a review of key investigation methods and results. We begin by describing the main assumptions in a QM/MM approximation. Next, we continue with an overview of the systems we are reviewing, and finally, we present for every chosen system a group of QM/MM results that we consider to be interesting applications that complement and expand the insight gained by experiments. For this work, we choose systems according to research directions in which our group has interests; thus, these examples are not exhaustive. Rather, they are representative of the knowledge that QM/MM methods can bring to the study of coinage metals interacting with biological and organic molecules.

2. QM-MM Method

The QM/MM method arrived on the scene in 1976, when Warshel and Levit proposed a method for mixing the advantages of quantum (QM) and classical (MM) approximations [32]. This method separates a smaller part of the complete system (CS), called the primary system (PS), and describes its electronic properties with quantum methods. The rest of the system, called the secondary system (SS), is simulated with classical approximations. Today, several QM/MM schemes exist that vary in how they describe PS and SS interactions [33–38]. Currently, the most used schemes are the subtractive and additive. Indeed, all the studies presented in this review use the additive scheme.

In the subtractive scheme, CS interactions are computed classically, then PS’s classical energy is subtracted and replaced by the quantum energy. This implies that, in this scheme, the classical part does not affect the electronic distribution. Therefore, the quantum part interacts classically with the rest of the system. In this scheme, the potential energy can be written as

(1)

where the subscript and superscript indicate the approach and the region, respectively. Since the complete system is described classically, this scheme requires a classical force-field that describes all the atoms which for metal-organic systems often is not available.

On the other hand, the additive scheme uses quantum approaches in PS and classical approaches in SS as independent systems, and adds an interaction potential between both systems. This is

(2)

where the last term is the interaction potential. This potential considers a hybrid electrostatic behavior defined by

(3)

with as PS’s electronic density, and the position and point charges of SS’s atoms, the atomic number of PS’s atoms, and where contains the non-electrostatic electronic interactions. In the Figure 1 a schematic comparison of additive and subtractive approach is presented.

Figure 1. (a) Example of a subtractive QM/MM scheme; (b) additive QM/MM scheme; and (c) zoom in to the QM/MM interface region of an additive coupling scheme.Reprinted with permission from E. Brunk and U. Rothlisberger Chemical reviews 115 (2015), pp. 6217–6263 [39]. Copyright 2015 American Chemical Society.

To apply QM/MM, we must define crucial aspects such as the border between the quantum and classical part [40,41] or the interaction potentials between PS and SS, which could require some ad hoc corrections [42,43]. These definitions depend on, for example, the existence of covalent bonds between atoms of both regions or the basis set’s accuracy in describing PS’s electronic distribution. A correct definition of these aspects could avoid typical errors as the ‘electron spill out’ problem [44].

There are several approaches to deal with the problem of cutting covalent bonds between PS’s and SS’s atoms [45–47]. The three most popular schemes are link-atoms, frozen orbitals, and pseudopotentials. The objective of these schemes is to mimic the effects of the cut covalent bonds over the electronic distribution and to cap the valence electrons to avoid artificial reactivity. The link-atoms scheme proposes adding artificial atoms to the PS (usually hydrogen atoms) [48,49]. The Frozen orbitals scheme fixes some orbitals that are not included in the self-consistent iterations [32,50]. In pseudopotential scheme, an additional term is added to the potential energy [51,52].

A complete explanation of these approaches is beyond this review’s scope, but are discussed in the references given and elsewhere [33,39,53,54].

3. Systems overview

Copper combined with organic ligands through Cu-N bonds can be excited by radiation to induce charge transfer states [55]. This property is limited by metal-to-ligand charge transfer (MLCT) state instabilities (described experimentally in other work [56–58]), but what remains unrevealed is these states’ complete relaxation mechanism. Here, we show efforts to characterize aspects influencing the relaxation process using QM/MM in the Copper(I) phenanthroline complex Cu(dmp) (Table 1, rows 1 and 2).

Table 1. Metal-organic systems studied with quantum and molecular mechanics (QM/MM) methods. The code, QM, and MM approximation are listed in the second, third, and fourth columns, respectively.

System	Code	QM	MM
[Cu(dmphen) ]	Amber 10 [17]	DFT-BLYP	MM [17,63]
	ASE-GPAW	DFT-BLYP + SCF	LJ-ElectroS
Atox1-Cu	CHARMM	SCC-DFTB	CHARMM27/TIP3P
Atox1-Cu-ATP7B	Siesta	DFT-PBE	Amber99/TIP3P
C Ag	CPMD	DFT-PBE-D3	AMBER
G -Ag -G	CP2K	DFT-PBE-D3	AMBER
C T +Ag	CP2K	DFT-PBE	parmbsc0/TIP3P
Ag @C -[Ag -C	CP2K	DFT-PBE-D3	AMBER
(C-Ag -C)	CP2K	DFT-PBE-D3	parmbsc0/TIP3P
Au (SR)	CPMD	DFT-PBE	FF99SB/TIP3P
Au (SR)	CPMD	DFT-PBE	FF99SB/TIP3P
Au ZW	DALTON	DFT-CAM-B3LYP	UFF
Au (SCH )	GAUSSIAN	BLYP	AMBER

As we touched upon briefly in the introduction, copper is essential in cellular processes such as respiration, iron oxidation, antioxidant defense, and connective tissue formation [59,60]. Investigating protein behavior and revealing molecular details can improve their understanding and guide how to deal with diseases associated with a failure in the biological mechanisms’ proper behavior. Copper in large biomolecules is present as Cu with typical Cu-S bonds. In this review, we present some details of the secretory path of Cu elucidated from a QM/MM point of view (Table 1, rows 3 and 4).

Relatedly, in the presence of DNA strands, silver ions in solution form stable complexes with interesting photoluminescent properties [4,5]. Both Ag and Ag clusters have been reported experimentally, but only recently in 2019 has a crystal-determined atomic structure been obtained [61]. For many years, therefore, QM/MM methods were an important source of insight of such systems’ geometric and electronic structure. The metal-organic bond is Ag-N, which can lead in the case of DNA-Ag clusters to various charges in addition to the expected Ag and neutral silver. Thus, here we also review the QM/MM studies focused on determining the atomic structure of Ag-homopolymers of cytosine and guanine and DNA-stabilized silver clusters (Table 1, rows 5 to 9).

Similarly, thiolate-protected gold was synthesized by chemical methods as early as the 90ʹs [6]. But only in 2007 was the first crystal of a thiolate-protected cluster (Au (SR) ) determined [62]. The crystal opened the way to the first ab initio study of their electronic structure, and then the super atom model was proposed [7]. As the model was determined by DFT ground-state calculations, all effects related to temperature and solvent were left out of the model. We review here a series of QM/MM studies on thiolate-protected gold clusters that address these effects (summarized in Table 1, rows 11 to 13).

Solvent effects not only modify these thiol-protected gold clusters optical properties but also affect bare silver clusters in noble gas properties. We included studies at the QM/MM level that show the medium effectively shifts peaks in the optical absorption spectrum.

4. Copper(I) phenanthroline complex

To begin, we look at a copper(I) phenanthroline complex [Cu(dmphen) ] , which has important applications in photoabsorption of ultraviolet-visible light that excites the system to an MLCT state. With this property, this complex is a promising candidate in dye-sensitized solar cells (DSSCs) [64]. However, this MLCT transition includes fast distortion and flattening of the ligands (Figure 2), which exposes the metal binding site to interactions with the solvent molecules and increases the MLCT state’s instability. This limits the use of the [Cu(dmphen) ] complex in practice.

Figure 2. Structures of ground (S ) and photo-excited (S ) states of [Cu(dmphen) ]. Blue and red vectors show the orthogonal direction of N-Cu-N planes. is the angle between these vectors and its reductions represents the flattening of ligands. Color code: Cu = pink, N = blue, C = gray, H = white. Reprinted figure with permission from [67].

The complex’s flattening dynamics can be reduced by ligand substitution, but a proper proposal of substitution requires a vast knowledge of the flattening mechanism, which is still under construction. Indeed, the changes in the conformation’s geometry and its behavior have been associated with different effects. On one hand, the distortion has been considered as a consequence of a Jahn-Teller (JT) instability [57,58]. On the other hand, Tahara et al. [56] proposed the existence of a shallow barrier of energy between both states.

Since the first studies [56–58], it was clear that the changes in the complex’s configuration caused by excitation in the short timescale is independent of the solvent. However, the complex’s lifetime presents a significant difference, depending on the solvent ( 1.6 ns in acetonitrile; 90 ns in dichloromethane). This dependence was first explained by considering that the solvent molecules approximate to the copper binding site once the excited complex is flattened, forming a solute-solvent exciplex, which is a stable state that appears in the excited state, although it does not exist in the ground state [65].

In this sense, Penfold et al. [66] studied the influence of dichloromethane and acetonitrile solvent in the flattening process. In this scenario, they considered the complex as the quantum part and the solvent as the classical one. From MD simulations, the authors assessed the time that water molecules linger in closest proximity to the complex’s copper atom. They found that the closest proximity of a solvent molecule remains in a range of 3.5 to 4 angstrom (Å), for a maximum time lapse of 200 femtoseconds (fs). Considering these results, the authors consider it infeasible to create an exciplex between the complex and solvent molecules.

Levi et al. [67] used QM/MM to present a complete analysis of the different aspects related with the flattening of [Cu(dmphen) ] in acetonitrile. They confirmed the absence of exciplex formation by considering the radial distribution of solvent atoms from the copper atom, obtaining a clear orientation of the acetonitrile molecule through the nitrogen atom–but not special changes in the concentration of solvent molecules close to the copper atom in the flattened configuration. This result agrees with Penfold et al.’s analysis [66].

Moreover, Levi et al.’s work presents an estimate of the minimum energy path between the excited orthogonal and flattened states using the Nudged Elastic Band (NEB) method. With this calculation, no energetic barrier exists between those states. This result increases the reasons to consider the JT effect to explain the flattening change.

Finally, they showed the relationship between the flattening and orientation of methyl groups of the complex, fixing the distance to the hydrogen atoms. They obtained a clear correlation of the methyl orientation and the flattening state, and a difference of the orientation in the presence or absence of solvent molecules. With these results, they suggest a possible way to avoid complex flattening by changing these residues in the complex, opening a potential avenue to improve possible applications.

5. Antioxidant 1 chaperone and copper transfer

Next, as mentioned, we consider that copper plays an important role in cellular processes such as respiration, iron oxidation, antioxidant defense, and connective tissue formation [59,60]. However, copper concentrations must be regulated carefully. Excess copper in cells could deregularize proteins’ oxidation and produce reactive oxygen species (ROS). The homeostasis process that avoids the excess of copper concentration in humans follows a known path: the copper is introduced in the cell by the hCrtl protein, then it is transported by the antioxidant-1 (Atox1) chaperone to ATP7A/B proteins in the Golgi apparatus that secretes the excess [60].

Perkal et al. researched copper’s coordination state in the Atox1 dimer chaperone [68], taking copper and 43 significant atoms of the chaperone as the quantum part and the rest of the protein and water molecules belonging to the classical part. Additionally, six hydrogen link atoms were added in the covalent bonds between the quantum and classical bonded atoms to mimic the electronic distribution in the quantum part. The authors support their results with Electron Paramagnetic Resonance (EPR) experiments.

In their simulations, Perkal et al. observed a variation between an open and closed conformation in the chaperone, characterized by a greater or lower distance between the three main motifs in the chaperone, respectively (Figure 3). Another change in the structure that they found in their MD sampling was related to copper-sulfur coordination. The closed conformation presented a correlation with the four coordinated states, while the open conformation showed a correlation with either two or three coordination states by having more flexibility and more water molecules close to the copper site. The authors associate the metal binding site’s hydration with a possible protonation of Cys15, which favors the open conformation’s stability. The open conformation is associated with the copper transfer from the hCtrl protein to Atox1 chaperone, and from this last one to the ATP7A/B protein in the Golgi apparatus.

Figure 3. Averaged distances (Å) between motifs of two Atox1 monomers correspondent to three (gray numbers) and four (black numbers) coordinated copper. Color code: brown = Cu, yellow = short -helix of Atox1 dimer, dark gray = large -helix of Atox1 dimer, cyan = sheet strands of Atox1 dimer. Reprinted figure with permission from [68].

The copper transfer mechanism from Atox1 to ATP7B protein also has been investigated using QM/MM. Rodriguez et. al [69] assessed the potential energy of all possible intermediates in the Cu-transfer reaction from the 2-coordinated Atox1-Cu(I) reactant to the 2-coordinated ATP7B-Cu(I) final product. The transition through each intermediate reaction was possible by adding a restrictive potential in terms of a reaction coordinate and equilibrating the system in each part of the reaction path. In this way, they proposed the most probable reaction mechanism. The proposed reaction is kinetically feasible because the maximum activation barrier is of 9.5 (kcal/mol). However, the resulting product is more energetic than the reactant, with an overall energy difference of 7.7 kcal/mol. Then, the authors suggest that thermal effects should contribute to the heterocomplex dissociation, such that the reaction is directional in vivo.

Details of Cu routes in biological processes–for example, the C(I)-C(II) reduction–remain unclear [60,70]. In silico analysis helps fill these gaps; likewise, the computational efforts we discussed here shed light upon characterizing the chaperone Atox1ʹs molecular structure with the copper atom and the Atox1 to ATP7B Cu transport.

6. DNA-stabilized silver clusters

Silver clusters stabilized by DNA complexes are nanobiosystems with high fluorescence, good photostability, low toxicity, low cost, and potential for novel biomedical and technological applications [4] with only one of these systems crystallized recently [61]. In this context, simulations and the QM/MM work reviewed here have played a pivotal role in the understanding the relation between its electronic properties, stabilization factors, and its geometry at atomistic level.

Espinosa et al. showed that it is essential to include the solvent in modeling the Ag-DNA atomic structure [71]. From simulations, they also proposed the existence of a novel DNA H-bond termed interplanar H-bond that participates in its stabilization. In that work, both the QM/MM simulations and the experiment focused on exactly the same system: a C Ag tetramer.

The tetramer complex can be described as the structure resulting from the Ag stabilization of two strands of DNA, each containing two cytosine nucleobases. Importantly, the simulated tetramer stabilizes with two interplanar H-bonds in solvent, while in vacuum with only one bond. When only one interplanar H-bond is present, the structure becomes wider and water would be able to interact with nucleobases in the central region. Thus, the compact two interplanar H-bond structure seemed a better candidate. Circular Dichroism simulations with time-dependent DFT (TDDFT) [72] and measurements indeed showed quite good agreement between the compact structure and the experiment. To avoid spurious transitions in metal-organic systems coming from generalized gradient approximation (GGA) exchange correlation functionals, the optical simulations were performed with Gaussian code and the long-range corrected functional CAM-B3LYP [73] on a geometry extracted from a QM relaxation of the structure surrounded with a shell of water.

Here, it is important to remember that there are at least two problems related to the simulation of optical properties with hybrid functionals. One is the largest computational time required compared to generalized gradient approximation (GGA) functionals, and the other is the possibility of having spurious transitions, depending on the amount of exact exchange added. A solution proposed for the former problem is in Makkonen’s work, using QM/MM MD trajectories obtained with PBE functionals [74] to average the simulated Circular Dichroism spectra [75]. Applied to the silver-stabilized guanine homopolymer G - Ag -G , it showed excellent agreement between a parallel left-handed-oriented tetramer and the experiment [76].

The hypothesis that silver-DNA complexes are indeed a mixture of previously studied Ag cations and clusters was used by Ramanazov et al. to determine a series of Ag @DNA and Ag @DNA violet-, green-, and red-emitting clusters [77]. The proposed structures were equilibrated in QM/MM MD simulations, and both the structure and charge were refined to match the experimental optical spectrum reported in the same study. The work proposed a partially oxidized cluster as the C T +Ag system’s structure.

In a similar way, Chen et al. started with a mixture of silver in single cationic form and aggregated in a cluster form stabilized inside DNA duplexes [78]. Several structures appeared equilibrated, and those whose chiroptical signal simulated via MD averaging matched the experiment was then studied in more detail [78] (Figure 4). An important result is the evidence that different Bader charges [79] can coexist in the cluster (Figure 5), depending on the exact coordination, and in agreement with previous proposal [77]. The partially oxidized cluster proposal appears in contrast to thiolate-protected gold or silver clusters, where the cluster’s metallic core only has a very small positive charge, which is indicative of a covalent S-Au [80].

Other QM/MM studies include the dissociation of the cytosine-Ag cations bonds [81] through free energy estimation via thermodynamic integration. The Free Energy landscape shows that the tetramer (C-Ag -C) has as global minimum: the transoid configuration with two stabilizing interplanar hydrogen bonds (H-bonds), in agreement with experiments and Espinosa et al.’s simulations. It shows an additional local minimum at longer cytosine-cytosine distances that may be relevant to understand experimental complex formation.

Figure 4. (a)–(c) Structure of the Ag @C -[Ag ] -C , Ag @C -[Ag ] -C ,and Ag @C -[Ag ] -C ; (d)–(f) the Ag core and their neighbor atoms in the Ag:DNAs, respectively. evidence. Reprinted figure with permission from [78]

Figure 5. (a)–(c) Bader charges of the Ag atoms in the Ag @C -[Ag ] -C , Ag @C -[Ag ] -C ,and Ag @C -[Ag ] -C .The top Ag atom was marked by 1 and the bottom Ag was marked by 6, 7,and 8, respectively. The labeling of the Ag core atoms can be seen in the inset of the figure. Reprinted figure with permission from [78]

7. Thiolate-protected gold clusters

A thiolate-protected gold cluster Au (SR) electronic superatomic structure appears as a closed shell with 8 electrons in S and P shells as 1S 1P . The unoccupied levels (also superatomic) are type D [7,82]. The first transition in the optical spectrum is the allowed P D transition. The ligands disposed in tetrahedric symmetry around the metallic core break the spherical symmetry and impose a splitting of the D degenerate orbitals. The first peak in optical spectrum is therefore separated in two peaks. Ligands SR in simulations were often chosen to be R = CH for convenience. Clusters also were placed in a vacuum with ground-state optimizations, followed by time-dependent TDDFT optical absorption simulations. Earlier work [82] noted an excellent agreement between the first simulated excitation with a GGA functional and experiment. It would have been expected, however, for the simulated gap to be underestimated with respect to an experiment taking into account GGA delocalization problems.

The QM/MM work from Rojas-Cervellera et al. [83] explains that the match between the absolute energy values of simulated and experimental absorption peaks comes from a cancellation of errors. QM/MM simulation of Au protected by gluthathione ligands solvated in water show a reduced overall atomic symmetry, which in turn expands the metallic core. The reduced symmetry splits even more the two D shells in the electronic structure, which then implies a reduction of the highest and lowest occupied molecular orbital (HOMO-LUMO) gap. Using both a realistic ligand and the solvent decreases the HOMO-LUMO gap, which then implies a reduction of the first excitation optical peak. Furthermore, the functional’s change to a more accurate hybrid functional (PBE0) [84] then drastically increases the HOMO-LUMO gap again, finally approaching it to the expected experimental energy range (Figure 6).

Figure 6. HOMO−LUMO gap of these systems with PBE and PBE0 exchange-correlation functionals. Green bars correspond to Au (SCH ) , whereas yellow bars refer to Au (GSH) . Reprinted with permission from [83]. Copyright 2015 American Chemical Society.

Thiol-protected Au (SR) was studied by Cheverier et al. [85] to elucidate how two different solvents affect the atomic and electronic structure. They found that in toluene, there is an increase in disorder that opens free space on the core surface. The surface core can then be reached for catalysis. Clusters with glutathione ligands, however–because of the glutathione–are much bulkier ligands. Independently on the solvent, they do not allow a solvent’s entrance to the core surface. They also found increased charge transfer from a glutathione ligand toward the core, which then could explain the observed optical activity enhancement.

A similar effect is found by Peric et al. simulating zwitterion functionalized gold clusters [86]. The ligands create an enhancement of luminescence via polarization. We want to note an interesting variation where only four ligands are included in QM, and the rest of the ligands and solvent are in MM. They propose a rigidification of the core appears when increasing the number of ligands. The ligands also remove the water solvent molecules from the surface and ligand region to the core’s outer region. Based on measured and simulated one- and two-photon absorption spectra, they propose using this ligand to increase its photoluminescence quantum yield.

Thus far, all the QM/MM regions we reviewed here have been chosen to include the metal on the QM part. There is, however, one work that proposed including the metal core in the MM region and the ligands in the QM region [87]. Applying the method to nuclear magnetic resonance (NMR) shows that an error with respect to the experiment is the same for QM/MM and a full QM.

8 Bare silver clusters in noble gas

As a final example, here we include simulations of the effect of a medium-like noble gas on the cluster’s optical spectrum. Earlier simulations [88] chose to place Ag and Ag in the QM region while placing He as explicit atoms in the MM region. The code used was Demon2k and the QM/MM coupling method was mechanical. The simulated and experimental spectrum of Ag agrees in fine detail. In the case of Ag , however, newer calculations with a more accurate model showed different effects [89].

A different choice of regions includes the metal and some noble gas atoms in the QM region while the rest of the solvent is present as an implicit solvent in the MM region [89] (Figure 7). Strictly speaking, it is rather a QM/classical method where the classical region is represented by an extended field. However, the combination of methods showed a high precision, and we think it is worth including here. The code used in the work is Gaussian. They obtained that from Ag ; and in clusters with higher number of atoms, there is a redshift in the absorption spectrum coming from the environment’s polarization effect. The higher the dielectric constant, the higher the shift that varies smoothly. In small clusters such as Ag8, there is a superposition of the classical redshift and the quantum effects, implying that some peaks are redshifted, others are blueshifted, and then others remain unchanged with respect to vacuum simulations.

Figure 7. Absorption spectrum of Ag in a neon matrix calculated with several dielectric constants (1.1–1.9) together with the spectrum calculated in the gas phase (green line) and the experimental spectrum in neon matrix 19 (black line). On the right, a scheme of the model showing Ag surrounded by 100 rare-gas (RG) atoms while the effects of the more distant atoms are modeled via a dielectric medium characterized by the dielectric function . Reprinted with permission from [89]. Copyright 2018 American Chemical Society.

9. Conclusion

The hybrid QM/MM method can complement experiments, revealing details of excited state evolution, charge transport, light absorption and emission, and the atomic structure in the absence of a crystal-determined structure. With that in mind, here we highlighted useful applications of QM/MM methods to study copper, silver, and gold atoms and clusters interacting with biological and organic molecules.

One particular aspect that we were able to show is how QM/MM simulations strengthen the hypothesis of the Jahn-Teller effect rather than forming a solvent-solute exciplex for explaining MLCT states’ short lifetime in [Cu(dmphen) ] . Moreover, we explored how QM/MM can be used to establish the relationship between the orientation of methyl groups and the flattening of the complex. Another useful application entails transporting a copper ion from one protein to another in human cells’ Cu secretion path. A fourth example shows studies of Atox1 copper states characterizing the flexibility in the dimer and redistributing relevant solvent molecules in the copper binding site.

QM/MM simulations also provide important structural information on silver-DNA complexes’ stabilizing factors. For example, a novel interplanar H-bond’s existence was proposed and then later confirmed via X-ray crystallography. Identifying clusters with mixed atomic charges, or oxidation, has been proposed in multiple QM/MM simulations.

Yet another fascinating use of QM/MM surrounds the fine detail that has been gained on how solvents and ligands affect the optical and geometrical properties of thiolate-protected gold clusters: both toluene and water increase the atomic geometrical disorder, which leads to a decrease of the HOMO-LUMO gap. The use of bulky ligands such as glutathione and zwitterion ligands can prevent the entrance of solvent molecules to the cluster surface and increase the metallic core’s symmetry. In the case of the zwitterion ligand, this process also correlates to a higher photoluminescence quantum yield.

Despite the great strides and bounty of successes that have resulted from using QM/MM, each of these methods are algorithms that still require large computational resources. They would benefit therefore from advances in faster quantum methods such as orbital-free DFT (OFDFT), potential energies from neural networks, and automating parameters related to the QM-MM cut or border. We expect that QM/MM methods will continue to boost supporting the exploration of novel hybrid organo-metallic materials and help us understand how these materials’ interplay gives rise to such diverse electronic properties.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Universidad de Antioquia and MinCiencias.