The mechanism underlying the functionalisation of cobalt nanoparticles by carboxylic acids: a first-principles computational study

The promise of biocompatible magnetic nanoparticles with high magnetic saturation in the implementation as drug carriers and hyperthermia agents has generated significant interest in functionalised cobalt nanoparticles. Carboxylic acid coatings on metallic nanoparticles have been shown as an attractive option owing to their respectable stability and biocompatibility. However, only limited information is available on the molecular mechanism leading to the formation of such protective coatings. In this study, ab initio molecular dynamics simulations have been used to unravel the functionalisation mechanism starting from a neutral cobalt cluster and valeric acid molecules. Three stages were detected in the coating process: (i) rapid initial adsorption of acid molecules, (ii) simultaneous adsorption of new molecules and dissociation of those already interacting with the cluster, and, finally, (iii) grouping of dissociated hydrogen atoms and subsequent desorption of acid molecules. The fate of the hydrogen atoms was probed through a combination of static and dynamic ab initio modelling approaches, which predicted H 2 generation with favourable energetics. A better understanding of the functionalisation and interaction mechanisms will aid the rational design of biocompatible cobalt nanoparticles for various applications.


Introduction
][3][4][5][6][7][8] Magnetic nanoparticle hyperthermia (mNPH) has been a promising player on the frontline of cancer treatment for about a decade, with many NP materials for this therapy being overshadowed by highly biocompatible iron oxides. 9,10However, despite much research into shape and composition, the advantage of the biological harmlessness of magnetite, maghemite, and many other metal oxides lost its magnitude under the ineffectiveness of insufficient magnetisation, which has fostered interest into highly magnetic metal NPs. 11,12Originating from the popularity of innocuous iron oxides, Fe NPs were a natural first choice and their smaller size and improved magnetic properties have suggested that iron-based magnetic carriers would be superior to the existing maghemite and magnetite NPs. 13 Combinations of Fe with other metals, such as FePt, [14][15][16] FeNi, 17,18 and FeCo [19][20][21] nanoalloys, also appeared attractive.However, the high magnetic moment of cobalt NPs gives them an intrinsic advantage over the aforementioned materials, which can be exploited specifically in mNPH, where a combination of magnetisation and small diameters is crucial for the desired performance.][24] Nevertheless, considerable challenges exist in adapting these metallic systems into pragmatic, functional biomaterials.The most important practical question regarding the design of metallic NPs for use in biomedicine is how to ensure a sufficiently long lifetime of NPs in the organism.This not only requires that the NPs do not disintegrate, aggregate, or lose their integrity otherwise, but also that they are not removed by the immune system.Additionally, for mNPH the NPs should maintain high levels of magnetisation throughout the duration of the therapy.These are significant impediments based on the interactions amongst NPs, as well as between the NPs and the biological environment, and it is suggested that they can be solved through NP functionalisation prior to administration. 25,263][34][35] Only a limited number of studies considered other types of ligands, such as polymers. 36,37Acid-encapsulated Co NPs have been tested for heating applications, including mNPH, 23,38 and as contrast agents for MRI, 39 and the results of these studies suggest that protected Co NPs show great potential as an alternative to iron oxides.
However, the multi-parameter nature of the stability of functionalised nanocomposites in respect to the physiological environment and phenomenological explanations at the atomic level are hardly accessible experimentally.The impact of the interplay between the NP core and ligands on the properties of coated NPs [40][41][42] makes it all the more important to fully understand these mechanisms and to be able to predict the structure of a proposed organic coating-magnetic core design.4][45] Gold has gathered by far the most interest for biomedical applications, owing to its inertness, radio-sensitizing capabilities, and easy functionalisation with thiols.Information collected on Au NPs both experimentally and theoretically ranges from the dependence of structural changes on the size of NPs, 46,47 through the mechanism of formation, [48][49][50] layout, 51,52 and influence of thiolate coatings on the membrane penetration, [53][54][55] to the irradiation of Au-SH nanocomposites and their performance in radiotherapy or drug delivery. 56Despite their promising qualities, Co NPs remain under-explored, with only a few studies on their reactivity and functionalisation.Our recent work has revealed that size and strain effects govern the alternations of geometrical shape with the increase in size of Co NPs. 57The follow-up study was focused on resolving contradictory experimental findings on the nature of the ligand-NP interactions in acid-functionalised Co NPs, assigning bridging bidentate as the preferential binding mode. 32,33,58However, virtually nothing is known regarding the functionalisation mechanism or reactivity of coated Co NPs beyond this basic information.The ability of the chosen ligand to function according to the biological requirements depends largely on the structure of the coating which is a complicated function of the surface density and type of the molecule in combination with material, size, and shape of the NP core, whereas the functionalisation mechanism can be a determining factor.
To explore the mechanism behind the functionalisation of Co NPs with carboxylic acids, we have performed ab initio molecular dynamics (AIMD) simulations of the process of coating with the molecules of valeric acid.The structure and stability of the formed nanocomposite were traced throughout the simulation with a particular focus on the changes in the adsorption rate and binding modes of adsorbed molecules.After the three stages of the functionalisation process were identified, the detected reactivity changes caused by the increased acid coverage and carboxyl hydrogen dissociation paved the way towards simulations of H 2 generation and desorption.From the calculated energy barriers, it appears that the formation of hydrogen gas can be expected within the process of carboxylic acid functionalisation of Co NPs, similar to the results of studies carried on thiol-protected Au NPs.

Models and methods
The characteristic length scales of functionalised NPs are of the order of a few nanometres with the number of atoms often beyond a few thousands.This puts constraints on the simulations, as the complexity of the system in terms of size by far exceeds the practical limits of precise ab initio methods which explicitly treat the electron dynamics of the system.With the interest of spanning the dynamics of the ligand molecules while also gaining electroniclevel information to monitor the energetic effects, ab initio molecular dynamics offered the best compromise in terms of system size and accuracy.
Inclusion of regions with increased activity, such as edges and vertices, through a cluster model increases the accuracy when predicting the dynamics of ligand molecules, compared to the limited movement possibilities on a flat surface.The size of the NP model was determined to mimic behaviour of experimentally realistic systems as closely as possible, while staying within the size ranges that can be simulated by the AIMD approach.For the purposes of functionalisation of Co NPs, a 57-atom cluster (approximately 1 nm in diameter) was shown to be at the lower limit of large nano-particulate behaviour regarding both electronic properties and carboxylic acid adsorption process. 58Therefore, in this work, 57-atom hcp Co cluster, with two (0001) and twelve (10% 11) facets, was utilised to represent the Co NPs, and valeric acid with its 5-C atom chain, CH 3 (CH 2 ) 3 COOH, was used as a ligand.
As a baseline for simulating the functionalisation mechanism, the reference case of a single particle in the presence of acid molecules was established.The previously optimised 57-atom Co cluster was placed in the simulation box together with 14 molecules of valeric acid (corresponding to a coverage of one molecule per facet if all molecules were to be adsorbed).The real density of the liquid was inaccessible to model by AIMD.AIMD simulations were performed using CP2K 59,60 with the PBE 61 exchange-correlation functional and the Goedecker, Teter, and Hutter (GTH) Gaussian basis set. 62The chosen system was simulated in the NVT ensemble during 10.00 ps of Born-Oppenheimer molecular dynamics.A Nose ´-Hoover thermostat 63 was set at the synthesis temperature of T = 500 K, 33,34,64 and a time step of 1.0 fs was used for the integration of the equations of motion.
Hydrogen dissociation and the formation of H 2 were studied by both density functional theory (DFT) and metadynamics simulations.Nudged Elastic Band (NEB) calculations 65 were performed using the Vienna Ab Initio Simulation Package (VASP) 66 within the generalized gradient approximation (GGA) of the PBE functional.The DFT-D3 method with Becke-Johnson damping was used to include the long-range dispersion interactions 67 and the k-point grid contained only the G point.In order to gauge the relevance of the size of the NP which was captured for the thiol-gold system, 68 H 2 generation on the 57-atom cluster was compared to the same process on a 13-atom counterpart with well-established gas-phase stability. 69n addition, a comparison in energy barriers of carboxyl and thiol ligands was carried out.To reduce the computational cost, smaller acetic acid and ethanethiol molecules were employed in these calculations.All intermediate and transition states were fully optimised and Hessian calculations were performed to verify the existence of one imaginary frequency for the transition states.CP2K was used instead for metadynamics trajectories, where Gaussian hills with a height of 0.01 Hartree and a width of 0.05 were deposited every 50 fs.VMD was used for post-processing. 70Results

Functionalisation mechanism
To investigate the functionalisation mechanism of Co NPs coated with carboxylic acids, a single 57-atom Co cluster was surrounded by 14 valeric acid molecules and AIMD relaxation was performed in a duration of 10.00 ps.Starting arrangement and optimised geometry after 10.00 ps are shown in the top panel of Fig. 1.The bottom panel of Fig. 1 contains the stimulation timescale with a few intermediate configurations capturing key events.
Initially, the cluster and the randomly oriented acid molecules were separated by a distance of about 8 Å. Acid molecules immediately started to move towards the cluster and were progressively captured on the surface via a monodentate interaction, which is characterised by a single bond between one Co atom and one O atom of the -COOH group.Adsorption of the first molecule occurred after just half a picosecond.The molecule adsorbed onto the vertex Co atom of the (10% 11) surface through the -COOH oxygen atom which was bound to the carboxyl hydrogen.Shortly after, a second adsorbate molecule bound to another (10% 11) surface vertex, but this time the bond was established between the Co atom and the freestanding oxygen of carboxyl functional group.At 1.00 ps, a total of five molecules had already interacted with the cluster, one of them forming a bond with the (0001) surface vertex atom.Interestingly, regardless of the initial orientation, the acid molecules showed no preference for adsorption by one or other of the oxygen atoms of the carboxyl group; even if the -OH part was initially further away than the corresponding freestanding O atom of the -COOH group, some molecules interacted through the former.Within the next picosecond, half of the molecules attached to the cluster, mostly to vertex Co atoms, but still no pattern was observed regarding the oxygen atom through which the monodentate adsorption took place.There were also instances where an interchange of the oxygen atom interacting with the cluster would occur within the acid molecule.The coating at this point could have been characterised as isotropic, i.e. adsorbed molecules were found at isolated positions and no aggregation took place.Just before the simulation time reached 2.00 ps, one carboxyl hydrogen dissociated from its parent acid and less than 0.15 ps later a second dissociation happened.Each of the dissociated molecules formed two Co-O bonds with two available cluster atoms, thus adsorbing in a bridging bidentate mode.The dissociated hydrogen atoms effortlessly moved across the cluster through available bridge and hollow sites.As the simulation proceeded, within the next two picoseconds two more acids lost their carboxyl hydrogen atoms, while two previously unbound molecules adsorbed forming single Co-O bonds with the cluster.Even though a portion of adsorbed acid molecules was still interacting through the oxygen atom that was connected to the carboxyl hydrogen, dissociation happened in all four occasions on acids that were binding through the freestanding oxygen atom of the -COOH group.Dissociated acids were observed to quickly form bridging bidentate modes of interaction, typically within 0.25 ps after dissociation.In one instance, intermediate chelate state, where both oxygen atoms bind the same Co atom, was observed.At this point, the movement of dissociated hydrogen atoms started to be more localised, with a pair from the two earliest dissociations pulsing around the same vertex Co atom, simply interchanging amongst each others positions.This caused the nearby acid molecule, which was actively interacting with the cluster for over two picoseconds, to desorb at about 4.00 ps, followed by the adsorption of another acid molecule at an available site farther away.Grouping of the other pair of dissociated hydrogen atoms was also captured just after 4.50 ps on the vertex on the opposite side of the cluster.At the end of the 5.00 ps, nine acid molecules were attached to the surface of the cluster and four had dissociated and were interacting in a bridging bidentate mode.The remaining five molecules were adsorbed undissociated in a monodentate binding mode.No significant disruption in the shape of the cluster was observed thus far, as indicated by minimal changes in the radius of gyration, r g , shown in Fig. 2.
Except for the constant repositioning of released hydrogen atoms, no further events were observed within the next half picosecond when an additional valeric acid molecule was adsorbed.With almost 40% of Co atoms on the surface of the cluster now directly interacting with oxygen atoms of adsorbed valeric acid molecules, shape of the cluster started to distort.Throughout the change, captured average values of the radius of gyration were around 3.96 Å, Fig. 2. Just after 6.00 ps, another dissociation occurred and it took almost a full picosecond for this latest dissociated carboxylic acid to establish a bridging bidentate interaction with the cluster because of the pulsating hydrogen atoms that were limiting the availably of surface Co atoms.This event was followed by the dissociation of yet another valeric acid at 7.00 ps, at which point the newly released H atom quickly paired with the hydrogen atom from the previous dissociation.By the time this acid formed a bridging bidentate mode with two Co atoms, the shape of the cluster had become icosahedral, and initial r g values were retrieved, Fig. 2. The next dissociation was observed at 7.70 ps, followed by subsequent grouping of hydrogen atoms and the desorption, just before 8.25 ps, of yet another valeric acid molecule, whose adsorption was observed in the early stages of simulation (at 0.45 ps).Last dissociation happened just after 8.75 ps of the simulation time, where the final picosecond was marked as the longest period without the occurrence of any events, before the simulation was terminated at 10.00 ps.
At the end of the 10.00 ps run, ten acid molecules were attached to the surface of the cluster, eight of which had dissociated and were interacting through the bridging bidentate mode.The remaining two molecules were adsorbed undissociated in the monodentate binding mode, one bonding to the vertex Co atom and the other to the edge Co atom.One of the O-Co bonds was established through the oxygen atom interacting directly with the carboxyl hydrogen.That specific valeric acid molecule underwent multiple attempts to form a permanent bond with the cluster (average O-Co bond of interacting acid molecules being 1.9 Å), however, for the majority of the simulation time it was at a 2.2-2.4Å distance from the cluster's surface, giving place to dissociated hydrogen atoms.The adsorbed acid molecules were distributed evenly across the whole cluster and both (0001) and (10% 11) surface areas, with the exception of the vertex sites where the dissociated hydrogen atoms resided.
To evaluate the stability of the adsorption system, potential energy during the 10.00 ps run was plotted in Fig. 3.The figure illustrates a significant decrease in the potential energy of the system during the early stage of the simulation up to 1.00 ps followed by a relatively stable period of about 0.50 ps.Another noticeable drop in energy occurred at approximately 2.00 ps also followed by a stagnant trend and a more subtle decline that extends to 4.25 ps when the first occurrence of a loss in stability can be seen.After about 0.25 ps, the negative slope of the potential energy curve is restored.These alternations in stability occurred throughout the remaining time of the simulation, whereas the on-going decrease in the potential energy faded away after 7.00 ps, seemingly reaching equilibrium state in the number of adsorbed carboxylic acid molecules.Two additional striking drops at 7.10 ps and 8.75 ps correspond to the stability gained through the shape modification, and a further adsorption event as a response to the hydrogengrouping triggered desorption destabilisation.
As the potential energy trend displays, the cluster functionalisation process can be divided into three distinct stages.Stage 1, in which a decrease in the potential energy is solely due to the effective adsorption of valeric acid molecules across the surface of the cluster at a relatively constant adsorption rate (first two picoseconds).In stage 2, the rate of adsorption is slowed down, and fluctuations in the potential energy of the system are caused by hydrogen dissociations and the establishment of bridging bidentate binding modes (between 2.00 ps and 4.50 ps).Finally, stage 3 comprises of alternating destabilisation and restabilisation periods characterised by desorption of valeric acid molecules near the locations of dissociated hydrogen atoms and counter-adsorption/dissociation of previously non-adsorbed/ non-dissociated acid molecules at distinct sites of the cluster.In this final stage, the adsorption rate has reached the plateau, and adsorption is only induced in the instances of hydrogen-caused acid detachment to regain the equilibrium state of the system.
A schematic representation of the AIMD-predicted reaction mechanism is shown in Fig. 4. Considering that the available literature on the binding of carboxylic ligands to monometallic and bimetallic NPs and surfaces reports the bridging bidentate as the preferred mode of interaction in most cases with the occasional dominance of the chelate biding, [71][72][73][74][75] it is likely that the proposed functionalisation mechanism, or certain aspects of it, are also applicable to those systems.The true mechanisms for the acid functionalisation of other metallic NPs are, however, subject to future research efforts.
][78] Changing the ionic strength of the medium might be a good strategy to obtain higher coating densities within the collapsedsite defects on the NP surface, as demonstrated on the example of Au NPs. 79he observed behaviour, together with the comprehension that such organic coatings are often synthesised with the loss of hydrogen, indicates that the path of hydrogen should also be investigated.Throughout the simulation time, as more acid molecules lost their carboxyl hydrogen atoms, desorption started to occur, thereby destabilising the system.Since it is very unlikely that hydrogen atoms would detach from the cluster on the timescale of the AIMD simulations, a combination of NEB and ab initio metadynamics was chosen to model the occurrence and energetics of this rare event.

Fate of dissociated hydrogen
Very little is known about the reaction pathways of hydrogen atoms dissociated during the conjugation of organic surfactants with metal NPs.Even in the case of well-investigated thiol-gold systems, there are only a few very recent theoretical studies of the hydrogen transfer mechanism to support an experimentally suggested release of H 2 in the functionalisation process. 50,80 previous AIMD work 81 has shown that interaction of two thiols with a four-atom gold cluster results in the production of H 2 and formation of substructures containing adatom-like gold atoms pulled away from the cluster surface with almost linear RS-Au-SR bonding patterns.84 Later, the interaction of small gold clusters with one to 38 atoms was examined in a DFT study, which confirmed the thermodynamics of the H 2 generation and subsequent staple motif formation.68 To the best of our knowledge, no information is available on cobalt-organic surfactant systems.To provide a greater understanding of the functionalisation mechanism of cobalt NPs by carboxylic acids, nudged elastic band algorithm was employed to investigate the dissociation, transfer, and recombination of hydrogen atoms associated with the carboxyl group, and the study was further extended to thiol ligands for comparison.To reduce the computational cost of expensive NEB calculations, 2C-atom chain ligands were used, i.e. acetic acid CH 3 COOH, and ethanethiol CH 3 CH 2 SH, instead of the initially adopted 5-C atom counterparts.The effects of coverage and finite temperature were tested through metadynamics simulations on the original system.
3.2.1 Single-molecule adsorption.Adsorption energies of single molecules of acetic acid and ethanethiol at various available sites of the 13-and 57-atom Co clusters can be found in Table 1.Adsorption is spontaneous for all investigated positions, with the strongest interaction on the 13-atom cluster leading to the adsorption energies of À1.59 eV for acetic acid and À1.81 eV for ethanethiol.The 57-atom cluster is characterised by two (0001) and twelve (10% 11) surface areas, and adsorption was probed in available positions on both surfaces.Adsorption on the facet site of the surfaces was shown to be only half as strong as adsorption on the vertex Co atoms, with adsorption energies of À0.57/À0.53 and À1.08/À1.07eV for facet and vertex adsorption of acetic acid on the (0001)/(10% 11) surface, respectively.Adsorption energies of ethanethiol were found to be À1.13/ 0.90 eV for the (0001)/(10% 11) facet adsorption sites, and À1.19 eV for adsorption on either of the cluster's vertices.
Next, the reaction pathways from the adsorption of acetic acid and ethanethiol on cobalt clusters to the formation of H-Co-OOCCH 3 and H-Co-SCH 2 CH 3 complexes were calculated to determine the rate-limiting step.The reaction pathways, together with optimised structures, are shown in Fig. 5 for acetic acid and Fig. 6 for ethanethiol, on both 13-and 57-atom clusters.
For acetic acid, the first transition state involves in both cases the transfer of the carboxyl hydrogen from acid molecule over a Co-Co bridge towards the nearest available three-fold adsorption site on the cluster.The interaction mode thus formed has no cobalt atoms shared between the dissociated hydrogen and bonding carboxyl oxygen atoms.Considering the limited size of the 13-atom cluster, the only possible transition is over the clusters edge, with the energy barrier of 0.70 eV.Extended surface areas of the 57-atom cluster allow transfer of hydrogen along the same facet during dissociation, which considerably reduces the energy barrier to 0.15 eV in the case of the (0001) surface, and to 0.13 eV in the case of the (10% 11) surface.Once the hydrogen is dissociated from the acid molecule, it is relatively free to move around the cluster over sites which are not saturated by Co-O interactions.This is shown by transition barriers of the subsequent hydrogen movement which were calculated to be 0.09 eV on the 13-atom cluster, and 0.20, and 0.18 eV on the (0001) and (10% 11) surfaces of the 57-atom cluster.Upon dissociation, the acid molecule proceeds to interact in a bridging bidentate adsorption mode.
The transition state of ethanethiol dissociation is characterised by a similar transfer of thiol hydrogen atom over a Co-Co bridge, with the rest of the molecule tilting in the opposite direction to form a second S-Co bond with the nearest cobalt atom.The difference compared to the acetic acid case is that the hydrogen atom adsorbs at the closest three-fold hollow adsorption site by sharing one of the cobalt atoms which interacts with the sulphur.This transition state lies 0.26 eV above the undissociated ethanethiol-13-atom cluster complex, and 0.44 eV above the undissociated ethanethiol-57-atom cluster complex.Energy barriers for the movement of dissociated hydrogen are still relatively low, i.e. 0.22 and 0.33 eV for the 13-and 57-atom cobalt cluster, respectively.
The dissociation transition state lies significantly below the energy of separated reactants and, thus, this reaction should proceed easily, as seen in the AIMD simulations, where spontaneous dissociation of valeric acid was observed on the 57-atom Co cluster.

Two-molecule adsorption
To develop an understanding of the formation of H 2 from dissociated hydrogen atoms, a second adsorbate molecule is required.Adsorption of the second adsorbate onto the site which is already coordinated, such that the two molecules both bind to the same atom of the cluster, was shown to be preferred in the case of thiol adsorption on the gold cluster, and it was accompanied by the rearrangement into a staple motif. 68The formation of such motifs, as well as the influence of the adsorption sites of the second molecule on the generation of H 2 , were probed on the 13-and 57-atom Co clusters.
The reaction mechanism for the adsorption and dissociation of the second acetic acid molecule at both uncoordinated and coordinated sites, and recombination of the pair of dissociated hydrogen atoms on the 13-atom Co cluster is shown in the left panel of Fig. 7.The dissociation of the second acid molecule is spontaneous and the movement of the dissociated hydrogen atom is unhindered, similar to the adsorption of a single acid molecule.The transition barrier is 0.18 eV when the second acid is adsorbed at the uncoordinated site on the cluster, and 0.03 eV when it is adsorbed at the coordination site of the first acid.However, desorption of H 2 requires for two hydrogen atoms to interact within a single vertex cobalt atom, which is followed by a slight pull-out of the interacting Co atom.This transition state requires overcoming a rather large energy barrier of 1.76 eV when two acids are adsorbed on remote sites of the cluster, but 1.13 eV if they are coordinated to the same Co atom.Adsorption of the second acetic acid onto the uncoordinated cobalt atom (E ads = À1.25 eV) is favoured over sharing the interaction site (E ads = À0.65 eV) by 0.60 eV in the undissociated state, and by approximately 0.35 eV after dissociation.Once the hydrogen molecule is formed, the cluster-acid complexes have similar stabilities.However, the reaction energy for the production of hydrogen gas starting from the adsorption of the second acid at the uncoordinated site is À0.33 eV, which almost triples to À0.89 eV when one cobalt atom is shared between the two adsorbed acids.The 57-atom cluster offers many distinct uncoordinated adsorption sites, but since the vertex atoms showed significantly favoured adsorption, only vertex positions will be considered for the adsorption of the second acid.There is no obvious preference between the (0001) and (10% 11) surfaces when it comes to the adsorption of the second acid in a position that is remote from the cluster site already interacting with previously adsorbed acid molecule, with the difference in the E ads of only 0.06 eV (E ads (0001) = À1.20 eV and E ads (10% 11) = À1.26eV).Optimisation of the two acids adsorbed to a common vertex atom of the (0001) surface has led to repositioning of one of the acids, resulting in their adsorption more remotely at the cluster.Adsorption on the uncoordinated portion of the cluster remains favoured by 0.60 eV over the shared Co atom site on the (10% 11) surface (E ads = À0.67 eV), as it was the case on the 13-atom cluster.
The subsequent reaction mechanism for the dissociation of the second acetic acid and H 2 generation at the (10% 11) surface of the 57-atom cluster is shown in the right panel of Fig. 7.In contrast to to the 13-atom cluster, the dissociation barrier for the second acetic acid is 0.80 eV when adsorbed at the uncoordinated site and 0.55 eV when sharing the adsorption site with the other acid.These energy barriers are considerably  higher than the 0.15 eV barrier found for the initially adsorbed molecule, although dissociation of multiple acids was indeed observed in the AIMD functionalisation run.The subsequent movement of dissociated hydrogen remained relatively unconstrained with energy costs of 0.10-0.25 eV.The desorption and formation of H 2 upon the adsorption of two acetic acid molecules at remote sites of the 57-atom cluster is accompanied by an energy barrier of 1.15 eV.The barrier is thus reduced by 0.61 eV compared to the same adsorption arrangement on the 13-atom cluster.Structural distortion of the desorption cobalt site is also significantly less pronounced.Where two acids are coordinated to share the interacting cobalt atom, desorption of hydrogen atoms is energetically more restricted, with a barrier of 1.35 eV, which is 0.22 eV higher than on the cluster with 13 Co atoms.Another pronounced difference between the two cluster sizes is the stability of formed acid-cluster complexes, with the shared coordination site complex on the 57-atom cluster being 0.60 eV less stable compared to its counterpart with two acids at distinct adsorption sites.The reaction energy of H 2 generation in reference to the acid-cluster complex with one dissociated and one undissociated acetic acid is À0.010 eV and it remains almost unchanged at À0.012 eV when the molecules share the adsorption site.
Results for the adsorption at the uncoordinated site of the (0001) surface are fairly similar to the (10% 11) surface and a detailed description has thus been omitted.The energy barrier for the dissociation of the second acid is calculated at 0.63 eV and that of the hydrogen movement at 0.18 eV.Transition state for H 2 desorption is 1.17 eV higher in energy compared to the preceding intermediate state.
Finally, a third acid molecule was adsorbed onto the 57-atom Co cluster.Based on the previous findings, only the interaction on a distinct (10% 11) facet of the cluster was considered, with the first two acids not sharing their adsorption sites.The sole adsorption of another acid in an undissociated form did not alter the barrier of H 2 generation and desorption, which was found to be 1.21 eV.However, upon dissociation, the energy barrier was reduced to 0.65 eV with a reaction energy for hydrogen recombination of À0.18 eV, Fig. 8.
Hydrogen gas generation was probed for ethanethiol as well, at both uncoordinated and coordinated adsorption sites, with the reaction mechanisms presented in Fig. 9.The transition state barrier for the dissociation of the second ethanethiol molecule on the 13-atom cluster is 0.20 eV when adsorbed at the uncoordinated site, and 0.36 eV in the case of a shared Co Fig. 7 Reaction pathways and accompanying optimised structures of the adsorption and dissociation of second acetic acid and subsequent recombination of hydrogen atoms into hydrogen gas on the 13-atom (left) and on the (10 % 11) surface of the 57-atom (right) Co cluster.Red pathway corresponds to the adsorption of second acid molecule at the uncoordinated site on the cluster, while black pathway corresponds to the adsorption at the coordination site of the initially adsorbed acid.R, IN(1,2,3), TS(1,2,3), and P stand for reactants, intermediate and transition states, and products, respectively.Grey, teal, red, and white spheres represent cobalt, carbon, oxygen, and hydrogen atoms, respectively.atom interaction between the two molecules, which is similar to the dissociation of the first ethanethiol molecule.Likewise, the dissociated hydrogen atom can move around the cluster with minimal effort, the movement barrier being around 0.20 eV regardless of the adsorption site.The energy barrier and reaction energy for production and desorption of H 2 are, however, much more sensitive to the adsorption site of the second ethanethiol molecule.For hydrogen gas to form and leave the cluster, a barrier of 1.49 eV needs to be overcome if the second adsorbate molecule adsorbs at a remote adsorption site.The energetic cost increases to 1.65 eV when the two adsorbate molecules interact within shared cobalt atom.However, a difference in the reaction energy to obtain the products owing to the initially adsorbed state of the two molecules has greater impact; adsorption within the same site lowers the reaction energy from the endothermic value of uncoordinated site adsorption of 0.32 eV to the slightly exothermic value of À0.03 eV.Thus, compared to the acetic acid, the adsorption of two molecules of ethanethiol on the 13-atom Co cluster is more favourable for H 2 formation when they share the adsorption site.Formation of a staple motif is not as pronounced as on the gold clusters, but Co atom positioned between the two sulphur atoms is slightly dislocated from the cluster and the distance to the central atom is lengthened by 5.8%.
Adsorption energies for the second ethanethiol molecule on the 57-atom Co cluster at the uncoordinated and at the site of the initially adsorbed ethanethiol are À1.97 and À2.04 eV, respectively.However, the energetic preference for interaction within the same Co atom is lost upon the dissociation of the second ethanethiol molecule, which is characterised by a 0.42 eV barrier when the two molecules are adsorbed separately, but a 0.58 eV barrier if they are adsorbed at the same adsorption site.
The transition of hydrogen from one available three-fold site of the cluster to another is energetically less consuming, 0.12-0.17eV.The transition state for the desorption of the dissociated hydrogen atoms in the form of a hydrogen molecule is accompanied by energy barriers of 1.36 eV for the two ethanethiol molecules at separate adsorption sites, and 1.52 eV for the staple-like arrangement of the pair of ethanethiol molecules.The barriers are reduced by 0.10 eV in average in comparison to the 13-atom cluster.Nevertheless, the reaction energy for the hydrogen recombination process is much more favourable when two adsorbates are not directly  interacting within the same cobalt atom (À0.36 eV) than when they share the interaction site (0.03 eV), which is in a sheer contrast to the 13-atom cluster thiol complexes.

H 2 generation metadynamics
A 5.00 ps metadynamics run was carried out to completely explore the free energy landscape, with the distance between the two dissociated -COOH hydrogen atoms implemented as a collective variable to capture the transition barriers for H 2 generation more realistically.The starting geometry was taken from the AIMD functionalisation run, prior to the first acid desorption event.The observed H-H distances and the resulting free energy landscape as a function of the H-H distance are shown in Fig. 10.
Initially, the two hydrogen atoms were bridging the same vertex Co atom at a mutual distance of 2.76 Å, where they were observed to interchange positions for a significant period of time throughout the AIMD run, as depicted on the top left of the Fig. 11.The same interchanging trend was continued in the majority of the first picosecond of the metadynamics simulation, with H atoms 'orbiting' around the vertex Co atom through hollow and bridge sites.Just before the simulation time reached 1.00 ps, the H-H distance approached that of the H 2 bond and the newly formed hydrogen molecule detached from the cluster, also shown in Fig. 11.The system existed in this state for about 0.30 ps, followed by a dissociation of the H 2 entity and readsorption of the H atoms at a site of the cluster which did not coincide with the detachment point.The hydrogen atoms briefly maintained a shared type of interaction within the same Co atom before exploring distant parts of the cluster away from each other.In the period between 3.50 and 4.00 ps, additional dissociations of two carboxylic acids were observed, inducing independent movement of the two H atoms of interest towards remote, newly dissociated H atoms to form new grouping sites.At about 4.50 ps, the two H atoms of one of the new H-H pairs formed the transition state previously observed for the two initial hydrogen atoms of interest and desorbed shortly after from the cluster in the form of H 2 .This event was taken as the termination point of the simulation, and it is depicted in the top right area of Fig. 11.
The bottom panel of Fig. 10 depicts the free energy landscape as spanned over the H-H distance collective variable.The numerous minima captured for distances above 2.00 Å correspond to the positioning of two dissociated hydrogen atoms in various hollow and bridge positions available around the cluster.With energies between À0.28 and À0.53 eV, different combinations of hydrogen atoms situated in hollow and/or bridge positions presented similar energetic favourability.The alternating environment of available sites (proximity of edges, vertices, and adsorbed acid molecules) also affects the energetics of obtained minima.The energy barriers for the transitions between these combinations range from as low as 0.03 eV up to 0.50 eV, which confirms the low energy cost of hydrogen movement over the surface of the cluster, as predicted by DFT.A minimum with energy of À0.27 eV at a 1.32 Å H-H distance corresponds to the structure where the two H atoms have left favourable multi-fold positions and have each formed a single bond with the same Co atom, which resembles the transition state structures obtained by DFT.However, in contrast to the DFT geometry, the interacting Co atom was not pulled out from the cluster structure during the metadynamics simulation.Finally, another minimum is observed for the H-H distance of 0.75 Å, which corresponds to the bond length of a diatomic H 2 molecule, with an energy of À0.39 eV.Hence, the transition from the simultaneous interaction of both H atoms within the same vertex of the cluster to the generation and detachment of the H 2 molecule is predicted to have an energy cost of around 0.27 eV.This is a notably lower energy barrier than that calculated by DFT, which gives an indication of the influence that the presence of adsorbed valeric acid molecules and the experimental temperature conditions have on the strain and electronic properties of the cluster.
Metadynamics investigation of H 2 formation within the thiolprotected four-atom gold cluster captured a rearrangement of gold atoms, followed by an energy barrier of 1.21 eV for the desorption of the hydrogen molecule. 81However, a model of this size has only a limited capability to predict the behaviour of NPs.DFT calculations at 0 K have already shown significant differences in the energetics as the cluster size grows and a transition barrier of H 2 generation and detachment on a 12-atom Au cluster was predicted at 0.67 eV, around half that of the four-atom one. 68his trend correlates with the reduction in the transition barrier obtained when going from the 13-to the 57-atom cobalt cluster predicted by our DFT calculations.The lower energy barrier of Au clusters compared to Co clusters emanates from the formation of thiolate staple motifs, which were not found to be as favourable and as geometrically pronounced for cobalt.Additionally, it is well-known that hydrogen binds more weakly to cobalt than to gold clusters, 85,86 offering a compromise in the energy costs despite the failure to form highly stable staple motifs.Overall, the 57-atom Co cluster used in this study is considered to offer a fairly good representation for the NP behaviour.The energy barrier for the detachment of H 2 obtained in the metadynamics simulations is found to be lower than what has been reported by DFT calculations for gold clusters, which can be explained by the differences in the binding energy of hydrogen and the increased number of adsorbate molecules present.
[89] In addition, most favourable adsorption sites are also casedependent, and hence the fate of dissociated hydrogen atoms captured here for the acid-functionalised Co NPs is not expected to be maintained across the transition metal row.[92]

Conclusions
The mechanism of functionalisation of cobalt NPs by carboxylic acids was probed with AIMD, which has shown that the process can be resolved in three stages: (i) stage of sole adsorption, (ii) stage of concurrent dissociation of adsorbed molecules and adsorption of non-interacting molecules, and (iii) stage of alternating dissociation, desorption, and counter-adsorption processes with the stagnant potential energy and extensive aggregation of dissociated hydrogen atoms.The latter has been shown to prevent any further adsorption of acid molecules and hinder the formation of a protective coating of higher density, despite 40% of surface atoms still being accessible.
The investigation of H 2 formation and detachment from the cluster's surface revealed energetically feasible barriers for the generation of hydrogen gas within the functionalisation process and reduction of the barrier is observed for larger metal clusters.Furthermore, the desorption energy of H 2 from Co clusters has similar trends as the well-studied thiol-gold systems when predicted by DFT, while metadynamics simulations, which incorporated synthesis temperature and multiple ligand molecules, showed a significant reduction in the energy barriers.The insights into carboxylic acid functionalisation provided by this study will help design cobalt NPs suitable for a range of biomedical applications.
the presented results, including how to access them, can be found in the Cardiff University data catalogue at http://doi.org/10.17035/d.2020.0121356877.

Fig. 1
Fig.1Timeline of the 10.00 ps ab initio molecular dynamics run of the [57-atom Co NP -14 valeric acid molecules] system with initial geometry given on the top left, the final geometry on the top right, and seven key events featured at the bottom in the order of their occurrence, i.e. first adsorption (0.60 ps), first dissociation (1.60 ps), formation of a bridging bidentate (2.50 ps), pairing of dissociated H atoms (3.20 ps), first desorption (4.40 ps), further grouping of dissociated H atoms (5.75 ps), cluster adopts icosahedral shape (6.80 ps).Grey, teal, red, and white spheres represent cobalt, carbon, oxygen, and hydrogen atoms, respectively.

Fig. 3
Fig. 3 Potential energy as a function of time during the 10.00 ps ab initio molecular dynamics run on the 57-atom Co NP -14 valeric acid molecules system.

Fig. 4
Fig.4Schematic representation of the AIMD-predicted mechanism for the acid-functionalisation of Co NPs.Colour scheme is as follows: steps shaded in blue are spontaneous and stabilising, steps shaded in red are destabilising, and steps shaded in purple are hindered by an energy barrier, but, once overcome, they lead to further stabilisation of the system.

Fig. 5
Fig. 5 Reaction pathways and accompanying optimised structures of single acetic acid adsorption, dissociation, and subsequent hydrogen movement on the 13-(left) and 57-atom Co cluster (right).Two mechanisms on the 57-atom cluster distinguish between the adsorption on the (0001) and (10 % 11) facets, with depicted structures corresponding to the (10 % 11) surface.R, IN(1,2), TS(1,2), and P stand for reactants, intermediate and transition states, and products, respectively.Grey, teal, red, and white spheres represent cobalt, carbon, oxygen, and hydrogen atoms, respectively.

Fig. 8
Fig. 8 Reaction pathway and accompanying optimised structures of the adsorption and dissociation of the third acetic acid and subsequent recombination of hydrogen atoms into hydrogen gas on the (10 % 11) surface of the 57-atom Co cluster.R, IN, TS, and P stand for reactants, intermediate and transition states, and products, respectively.Grey, teal, red, and white spheres represent cobalt, carbon, oxygen, and hydrogen atoms, respectively.

Fig. 9
Fig.9Reaction pathways and accompanying optimised structures of the adsorption and dissociation of second ethanethiol and subsequent recombination of hydrogen atoms into hydrogen gas on the 13-atom (left) and on the (10 % 11) surface of the 57-atom (right) Co cluster.Yellow pathway corresponds to the adsorption of second ethanethiol molecule at the uncoordinated site on the cluster, while black pathway corresponds to the adsorption at the coordination site of the initially adsorbed ethanethiol.R, IN(1,2,3), TS(1,2,3), and P stand for reactants, intermediate and transition states, and products, respectively.Grey, teal, yellow, and white spheres represent cobalt, carbon, sulphur, and hydrogen atoms, respectively.

Fig. 10
Fig.10Top: H-H distance progression during the 5.00 ps metadynamics run for the system comprising of 57-atom Co NP and 14 valeric acid molecules; Bottom: Resulting free energy landscape as a function of the H-H distance as a collective variable.

Table 1
Adsorption energies, E ads , and Co-O distances, d Co-O , after the adsorption of acetic acid and ethanethiol on 13-and 57-atom Co clusters as calculated by DFT