The effect of ligands on the size distribution of copper nanoclusters: insights from molecular dynamics simulations

Controlling the size distribution in the nucleation of copper particles is crucial for achieving nanocrystals with desired physical and chemical properties. However, their synthesis involves a complex system of solvents, ligands, and copper precursors with intertwining effects on the size of the nanoclusters. We combine molecular dynamics simulations and DFT calculations to provide insight into the nucleation mechanism in the presence of a triphenylphosphite ligand. We identify the crucial role of the strength of the metal-phosphine bond in inhibiting the cluster's growth. We demonstrate computationally several practical routes to fine-tune the bond strength by modifying the side groups of the additive. Our work provides molecular insight into the complex nucleation process of protected copper nanocrystals, which can assist in controlling their size distribution and, eventually, their morphology.


Introduction
Metal nanoclusters (NCs) exhibit different physical and chemical properties as a function of their crystal morphology and size distribution. 11][12][13][14][15] In recent decades, significant advancements have been achieved primarily in the realm of noble metals such as Au and Ag by reducing a metal precursor in the presence of organic ligands. 9,16,179][20][21] Therefore, understanding the nucleation and growth of cooper nanocrystals in different environments allows for finer control over their morphology.This allows a monodisperse starting material for nanoparticle production or to effectively develop ligandprotected copper nanoclusters. 18,21Also, insights into nanocluster aggregation are necessary for the conventional colloidal synthesis of inorganic nanoparticles. 18,22nocluster synthesis usually involves adding ligands to precursors and solvents.Ligands can affect the size distribution and final shape of the NCs. 19Therefore, new ligands are experimentally investigated by adding them directly during NCs' synthesis. 23Ligands that form metal-phosphine complexes allow the bonding of the ligand group to the metal and thus alter the nucleation process in solution. 24However, the landscape of possible ligands is too vast to be explored only by a trial-and-error experimental approach. 25D simulations provide a microscopic view into the nucleation process as it evolves with time in atomistic detail, 26 which can, in principle, aid the design of ligands that lead to a desired size distribution. 25Realistic MD simulations of crystallization and nucleation from solution are challenging due to the multi-component system's complexity and its long-time scales. 27[46][47] In many experimental systems, ligands play a crucial role in controlling metal nanoclusters' stability and size. 48,49A recent experimental observation showed the significant effect of ligands on the resulting size distribution in the synthesis of copper nanoparticles. 22In that study, copper-ligand systems were investigated to control the crystallization and growth of copper nanoparticles.The synthesis procedure included mixing a precursor, solvent, and additive.One of the experimental procedures involves dissolving 0.5 mmol of copper nitrate in 8 ml Octadecene and 1 gr Hexadecylamine at 250 − 260 • C. Nucleation of copper in Octadecene showed large nanoparticles with an average diameter of 200 nm.However, when 1.6 ml of the ligand Triphenylphosphite (TPOP) was added to the solution, smaller nanoparticles were observed with an average diameter of 10 nm. 22 this paper, we use MD simulations and DFT calculations of this experimental model system to describe the ligand-protected copper nanocluster formation mechanism and the ligand influence on the resulting size distribution.We first provide the computational details and force field validation, followed by the nucleation simulation results.
We performed DFT calculations using the revPBE0 functional, 50 with the D3 dispersion correction and Becke-Johnson Damping. 51We employed the 6-31G(d,p) basis 52 and used Qchem 6.0 53 to perform a geometry optimization and find the equilibrium bond length between a Copper (Cu) atom and a Phosphorus (P) atom of the TPOP molecule.We applied a basis set superposition error (BSSE) correction 54 on the optimized geometry of the complex to find the bond strength.
We performed MD simulations using the Large-scale Atomic/Molecular Massively Parallel Simulator, LAMMPS (30 Jul 2021) 55 with the Nosé-Hoover thermostat and barostat 56 with a damping factor of 100 fs and 1000 fs for the temperature and pressure, respectively.
The simulations were conducted with a timestep of 1 f s.Lennard-Jones (LJ) parameters for non-bonding interactions between different atoms or pseudo-atoms were calculated using traditional Lorentz-Berthelot combination rules. 57We employed the united-atom TraPPE force field (FF) for Octadecene 58 with four types of pseudo-atoms: CH 3 (sp 3 ), CH 2 (sp 2 ), CH(sp 2 ) and an additional CH 2 (sp 2 ) type for the vinyl group.We tested this FF by calculating different properties in unbiased MD simulations at room temperature.We calculated the density in a simulation of a liquid with 1000 molecules and the bond lengths and angles in a simulation of a single molecule in a vacuum.The results showed a good comparison to experimental results (Table 1).GAFF2 FF parameters for TPOP were calculated using Pysimm for an all-atom description.We then replaced all carbons and hydrogen atoms with CH pseudo-atoms to reduce the computational cost.For these pseudo-atoms, we chose LJ parameters for the phenol CH group in the TraPPE FF.We tested our model using unbiased simulations of TPOP, with 1000 molecules in a liquid when calculating the density and a single molecule in a vacuum when calculating optimal bond lengths (Table 2).We used the embedded atom model (EAM) FF to describe the Cu-Cu interactions.To test this model, we evaluated the melting temperature of Cu using Metadynamics 63 simulations performed with PLUMED. 64We used the ENVIRONMENTSIMILARITY 65 collective variable with a bias pace of 500 timesteps, Gaussian height of 0.5 k B T and width of 3, and bias factor of 50.We retrieved the melting temperature of bulk Cu by performing simulations at different temperatures and calculating the free energy difference between the solid and liquid states (Figure 1a), showing that the EAM model provides a melting temperature of 1325K, close to the experimental value of 1357K. 66Interactions of Cu with other atoms were described with LJ interactions, found using Lorentz-Berthelot combination rules.However, when applying these rules, LJ parameters for Cu-Cu interactions are required.We used parameters as suggested by Lv et al. 67 We emphasize that these LJ parameters were used only for calculating Cu-solution interaction.Nonetheless, we verified that these LJ parameters adequately describe Cu by calculating the radial distribution function of bulk Cu (Figure 1b).We compared it to the embedded atom model (EAM) and calculations by Neogi et al.

Nucleation simulations
To test our hypothesis that chemical bonds between TPOP and Cu are responsible for changing the size distribution of the nanoparticles, we performed MD simulations of nucleation at the relevant experimental temperatures.In our simulations, which used relatively high temperatures, nucleation happens spontaneously, and Metadynamics simulations are not required.
The first simulation was performed in the isothermal-isobaric (NPT) ensemble at ambient conditions, including a premixed copper in solution without additive.The P-Cu bond strength was set to 15.3 kcal mol −1 , based on preliminary DFT calculation of the complex (see more details below).A snapshot from the simulations is shown in (Figure 2a), which exhibits nucleation to a single nanocluster.However, in the presence of the additive TPOP, significant multiple smaller nanoclusters are observed, as shown in (Figure 2b).Each copper nanocluster is surrounded by several TPOP molecules, as presented in (Figure 2c).The P atoms bond to the copper atoms, while the residues of the TPOP are oriented toward the outer side of the nanocluster.Thus, the TPOP molecules inhibit the growth of the copper nanocluster and allow a smaller copper cluster size distribution.The findings were corroborated by three independent trajectories.We tested the effect of the bond strength between the P and Cu atoms on the number of copper nanoclusters formed and their size, as shown in Figure 3.As the P-Cu bond strength decreases, the average number of atoms in the Cu nanoclusters increases, and the overall number of copper clusters decreases.Thus, altering the P-Cu bond strength can provide fine-tuning for the size distribution of the produced copper clusters.For example, modifying the ligand by adding an electrons withdrawing group can weaken the P-Cu bond.The weaker P-Cu bond allows copper clusters to grow and increase in size by merging smaller clusters.In addition, we tested the effect of temperature on the number of nanoclusters during the nucleation simulation.We performed simulations at temperatures of 190, 210, 230 and 250 • C. As the temperature decreases, the number of nanoclusters after 400ns slightly increases (Figure S1).Lower temperature decreases the diffusion of Cu atoms, thus hindering the formation of larger clusters, but the effect is less profound than the P-Cu bond strength.Using DFT calculations, we tested some practical approaches to increase the P-Cu bond strength experimentally.We looked at the effect of swapping two hydrogen atoms with electron-donating (-NH2) groups (Figure 4b).We also tested how swapping the phenyl residue (Figure 4c) with a smaller (methyl) residue affects the bond strength.We optimized the structure of the complexes with the suggested ligands and performed harmonic vibrational frequency calculations to confirm that they are local minima with no imaginary frequencies.Figure 4 shows the highest occupied molecular orbitals (HOMO) and the lowest unoccupied molecular orbitals (LUMO) for the original and altered complexes.
Swapping only hydrogen atoms with an NH2 group (Figure S2a-c) only slightly increases the bond strength (15.6kcal mol −1 ).Yet, when swapping two hydrogen atoms on differ-ent phenyl residues, as shown in Figure 5b, the bond strength increases to 19.8kcal mol −1 .
Swapping three hydrogen atoms with NH2 groups on different phenyl rings (Figure S2d-f) achieves a similar bond strength (20.5kcal mol −1 ).In the latter case, the molecular conformation of the complex is also very different relative to the original ligand (Figure 4a,b).
When replacing one phenyl ring with a methyl group, the bond between the copper and the altered TPOP increases even further to 28.8kcal mol −1 ,and the conformation of the complex changes (Figure 4a,c).Changing the phenyl ring with a smaller residue will also affect the evaporation of the ligand during the high-temperature synthesis process.Thus, only one phenyl ring was replaced to avoid significant changes in the physical properties.The suggested examples of modifying the ligand provide practical ways of controlling the complex binding energy and, as a result, the nanoparticle size distribution through additive design guided by MD simulations.

Conclusions
To summarize, we investigated the nucleation of ligand-protected copper nanoclusters in order to design ways to control their size distribution.We showed that adding a TPOP additive resulted in smaller NCs relative to a system without it.We used MD simulation to follow the nucleation process and understand the mechanism leading to control over NCs size.The metal-phosphine bond of copper atoms and ligand blocks the coalescence and growth of copper clusters, thus enabling smaller NCs.We showed that decreasing the P-Cu bond strength resulted in larger copper nanoclusters.P-Cu bond strength equal to or below 6kcal mol −1 resulted in a single nanocluster, similar to the nucleation process without the additive.We used DFT calculations to examine the effects of modifying the side group of TPOP on the P-Cu bond strength.Altering the ligand with such groups can allow finetuning of the copper nanocluster size, which is interesting to test experimentally.Our results show the potential of MD simulations in aiding the design of improved ligands for controlling the size distribution of nanoclusters with desired properties.

Figure 1 :
Figure 1: (a) Melting temperature calculation using EAM 69 for Cu-Cu interactions model.The orange line indicates the obtained melting temperature calculated at the intersection with ∆G=0 (b) RDF calculation for Cu-Solution interactions using LJ 70 and comparison to simulation reference.68

Figure 2 :
Figure 2: MD simulation of copper nucleation: (a) Copper nanoclusters without TPOP (for clarity, only copper atoms are shown), (b) Copper nanoclusters with TPOP (for clarity, only copper atoms are shown), and (c) snapshot of copper nanocluster and TPOP molecules (for clarity only Cu, P, and O atoms are shown in burgundy, orange, and red).

Figure 3 :
Figure 3: (a) Number of observed copper nanoclusters as a function of time for various P-Cu bond strengths, and (b) Mean number of atoms in a nanocluster after 400 ns of MD simulation, as a function of bond strength.Colors correspond to the same bond strength as in the legend of panel (a).

Figure 4 :
Figure 4: DFT calculation of TPOP additive with Cu atom, optimized geometry, HOMO and LUMO (a) original molecule, (b) modified TPOP with electrons donating group (-NH2), and (c) modified TPOP with methyl group.Atom colors of H, C, O, P, Cu, and N are white, gray, red, orange, gold, and blue, respectively.

Table 1 :
Comparison between experimental and simulated values of density and bond length for Octadecene 60experimental values obtained from ref59b experimental values obtained from ref60

Table 2 :
Comparison between experimental and simulated values of density and bond length for TPOP 62experimental values obtained from ref61b experimental values obtained from ref62