Theory and simulations of light-induced self-assembly in colloids with quantum chemistry derived empirical potentials

Abstract

Light-induced self-assembly (LISA) is a non-invasive method for tuning material properties. Photoresponsive ligands coated on the surfaces of nanoparticles are often used to achieve LISA. We report simulation studies for a photoresponsive ligand, azobenzene dithiol (ADT), which switches from a trans-to-cis configuration on exposure to ultraviolet light, allowing self-assembly in ADT-coated gold nanoparticles (NPs). This is attributed to a higher dipole moment of cis-ADT over trans-ADT which leads to a dipole–dipole attraction facilitating self-assembly. Singh and Jha [Comput. Theor. Chem., 2021, 1206, 113492] used quantum-chemistry calculations to quantify the interaction energy of a pair of ADT ligands in their cis and trans conformations. The interaction energy between ligands was fit to a potential energy function of the Lennard–Jones (LJ) form having distinct exponents for attractive and repulsive contributions. Using this generalized equation for the ligand–ligand interaction energy, we calculated the total effective interaction energy between a pair of cis as well as trans ADT-coated NPs. Specifically, we calculated the effective interaction energies between cis/trans-NPs using discrete as well as continuous approaches. Given the limitations of experiments in probing individual ligand conformations, we also studied the effect of varying the functional ligand length on the interaction energy between NPs and identified the optimal functional ligand length to capture the steric and conformational effects. Finally, using the effective interaction energy, we obtained a generalized potential energy function, which was applied in Langevin dynamics simulations to capture self-assembly in NPs.

---

1 Introduction

Self-assembly in colloidal systems can be controlled dynamically through various stimuli^1–5 including light,^6–9 temperature,^10,11 and pH¹² or electric and magnetic fields.^13,14 Light-responsive materials are particularly desirable due to several advantages of using light as a stimulus.^15–18 Light-induced self-assembly (LISA) leverages photoresponsive molecules, which respond to light by undergoing chemical or structural changes, to precisely control inter-particle interactions.^19,20 This enables a dynamic and reversible assembly with applications in materials, optics, and bio-medicine.^21–23 The advantages of using light as a stimulus are its noninvasive nature, high spatiotemporal resolution, applicability to a closed system with high specificity, and the flexibility to control the intensity and wavelength for different components in a system.

Experimental and computational studies have been carried out to understand LISA in colloidal particle systems focusing on various aspects like diffusion, nucleation, growth-control, and phase transitions.²⁴ The choice of the photoresponsive ligand plays a key role in determining inter-particle interactions to capture attraction, repulsion, or specific binding.^18,23 Some of the commonly used ligand groups to functionalize nanoparticles are azobenzene and spiropyran that undergo cis–trans isomerization reactions, or ligands like anthracene and coumarin which undergo photodimerization reactions.¹⁸ LISA in nanoparticles coated with azobenzene and its derivatives has been studied via experiments and theory.^25–30 Azobenzene molecule is a tunable ligand through modifications of functional groups attached to azobenzene.^9,31 The light-responsive nature of azobenzene ligands has also been utilized to drive the dissipative assembly in colloidal systems which find applications as self-erasing inks, tunable catalyts, and delivery vehicles.^32–34

Klajn et al.⁶ demonstrated LISA in gold nanoparticles coated with 4,4′-bis(11-mercaptoundecanoxy) azobenzene dithiol (ADT), which can switch between cis and trans configurations (Fig. 1A). We simply refer to them as NPs in this work. They studied self-assembly in these ADT-coated NPs by varying solvent composition (e.g., methanol, toluene, etc.) and the ligand coating density. Due to the hydrophobic nature of the trans-ADT ligand, the trans-NPs remain unaggregated in toluene since it is non-polar. However, the cis-ADT ligand develops a dipole moment, and in a polar solvent like methanol it aids in aggregation.³⁵ Therefore, the cis-NPs formed different structures like dispersed states, amorphous aggregates, reversible crystals, permanent crystals and supraspheres due to an interplay of temperature, solvent, and dipole moment.⁶ Due to a higher ligand density, irreversible permanent structures were also formed via covalent cross-linking of cis-NPs. However, at an optimum number density of ligands (∼20), the aggregation was reversible, and driven by dipole–dipole interactions between the azo-groups of the ligands. In trans-NPs, thermal energy dominates over dipole–dipole interaction and the particles remain dispersed. However, for cis-NPs, the dipole–dipole interaction can overcome thermal energy, thereby leading to aggregation of particles. Singh and Jha quantified this dipole–dipole interaction energy between a pair of ADT ligands in their cis and trans conformations.³⁶ They calculated the energy as a function of the separation distance between the center of mass of the azo-groups, by employing quantum-chemistry calculations performed using density functional theory (DFT). This interaction energy was fit to a generalized equation of the Lennard–Jones (LJ) form, but with different exponents for the attraction and repulsion in both cis and trans isomers.

Fig. 1 (A) Schematic representations of the chemical structures of the cis and trans isomers of ADT. The azo-group (blue) consists of two benzene rings connected via an N=N bond, and the alkyl-thiol chains (orange) consist of a 11-carbon chain with a thiolated (–SH) end. (B) Schematic representations of the cis-to-trans isomerization of a ligand-coated NP. The ligand-coated NP consists of 3 components: (1) NP core (yellow sphere), (2) capping agents DDA/DDAB (green lines on the NP surface), and, (3) the ADT ligand (blue azo-group and orange alkyl-thiol chains). Visible light stabilizes trans-NP forming a disorderd state, while UV light stabilizes the cis-NP leading to self-assembly.

Specialized models are required to study non-equilibrium thermodynamics of systems which undergo LISA. A scaling theory and a kinetic Monte Carlo approach has also been used to study LISA in a system of ADT-coated NPs.^37,38 The particles were modeled as spheres interacting via the LJ potential in the ‘ON’ state (UV light) and via the Weeks–Chandler–Andersen (WCA) potential in the ‘OFF’ state (visible light). This model was able to capture the general characteristics of the system with respect to aggregation dynamics.³⁷ Further, the effect of cyclic energy delivery was systematically studied showing that the periodicity of light exposure is more relevant for controlling self-assembly than a continued exposure of the system to light.³⁸ These previous studies described the aggregation kinetics of NPs under LISA. However, the generic models used previously are not suitable to explain the phase-equilibrium of the system, where interactions between the ligand-coated NPs need to be accounted for. For example, the range of interaction affects the phase-equilibrium, even for a simple square-well potential.³⁹ Therefore, a more accurate model that describes a pair of interacting NPs is needed for accounting inter-particle interactions.

The self-assembly of ligand-grafted particles is commonly studied using molecular dynamics (MD) simulations, considering ligands as coarse-grained beads on the NP surface and addressing the overall interactions between these beads.^40–42 This approach however is computationally costly and is not suitable for studies of larger systems. However, gaining inspiration from this method, we derived a generalized equation for the interaction energy for the ligand-coated NP system, incorporating the sum of individual ligand–ligand interactions contributing to the energy between a pair of ligand-coated NPs. The potential energy equation derived from DFT was applied to calculate the individual ligand–ligand interactions.³⁶ The generalized equation for the interaction energy was applied in Langevin dynamics simulations to briefly study self-assembly in ligand-coated NPs. The radial distribution function, average cluster size, and radius of gyration of self-assembled morphologies were analyzed. Several characteristics of the system like the minimum distance between aggregated cis-NPs, the type of clusters formed, and the soft repulsion in trans-NPs, are captured by the model.

2 Theoretical model and simulation methods

2.1 Description of ligand-coated NPs

We modeled the ligand-coated NPs based on the experimental system described by Klajn et al.⁶ Specifically, the NPs are considered to be of 5.6 nm diameter and coated with ADT ligands as well as a capping agent (dodecylamine, DDA) and a surfactant (didodecyldimethylammonium bromide, DDAB) for stability. Fig. 1B shows a schematic representation of a cis-NP and a trans-NP, and the reversible change in the conformation of the ligand-coated NP on exposure to UV and visible light. Each ligand-coated NP consists of three components: (1) NP core, (2) DDA and DDAB molecules attached to the NP surface, and (3) ADT ligands consisting of the azo-group and two alkyl-thiol chains. The azo-group of the ADT molecule consists of two benzene rings connected via an N=N double bond at the center. In the presence of UV light, ADT adopts a bent structure termed a cis conformation (Fig. 1A).^43,44 Upon exposure to visible light, ADT undergoes a cis-to-trans isomerization. The cis and trans configurations have a different spatial morphology. The cis-ADT molecule has a bent structure with its apex at the N=N bond. However, trans-ADT has a linear configuration. The photoisomerization mechanism of azobenzene has been studied rigorously and was found to occur through an ‘inversion-assisted torsion’ mechanism.^45–47 One thiol chain of ADT is attached to the NP surface, whereas the other is freely suspended. This free thiol end does not attach back to the same NP due to steric constraints around the N=N double bond.⁴⁸

In Fig. 2, we show the progressive simplification of a cis-NP and a trans-NP to the theoretical forms used in our model. The DDA and DDAB are neutral non-interacting groups and do not play a role in the attraction between the particles. However, they contribute to the excluded volume of the ligand-coated NP, and therefore are modeled as a non-penetrable shell around the NP core. The length of one thiol chain in ADT (11-carbon chain) is approximately equal to the length of one DDA/DDAB chain (12-carbon chain). Therefore, the azo-group of the ADT is considered to be present on the surface of the DDA/DDAB shell in our model. The azo-group in the cis-NP is exposed to the bulk solution making it accessible to the azo-groups of other ligands on other particles. The thiol chain attached to the NP surface is incorporated into the DDA/DDAB shell. The free thiol chain in cis-NP is neglected as it remains very close to the surface of the DDA/DDAB shell (labeled 2 in Fig. 2A). Based on these considerations, a conceivable model for cis-NPs consists of a hard-core sphere covered with a uniform distribution of points that represent the azo-groups (labeled 3 in Fig. 2A).

Fig. 2 Progressive simplification of (A) a cis-NP, and (B) a trans-NP, from the structural form to the theoretical form used in our model. (A) and (B) – (1) Representations of the cis and trans NPs consisting of all three components – the NP core, the DDA/DDAB capping agents, and the ADT ligands. (A) and (B) – (2) The DDA/DDAB is modeled as a non-penetrable hard layer (gray) around the NP. The azo-group of the ligand is present on the surface of the DDA/DDAB shell. The thiol chains at the free ends of the ligands are close to the DDA/DDAB surface for the cis-NP, but is extended away from the surface in the trans-NP. (A) and (B) – (3) The thiol chains attached to the NP surface are included in the gray layer. The free thiol chain close to the gray surface in the cis-NP is neglected. The free thiol chain extended away from the gray surface in the trans-NP is modeled as a soft-repulsive layer (orange). If the orange layers of two trans-NPs overlap, they experience a soft repulsion.

For the trans-NP, the free thiol chain extends out into the bulk solution, placing constraints on the closest approach of two interacting particles (labeled 2 in Fig. 2B). This restricts the azo-group of one particle from being in proximity to a similar group on another particle. The additional length of the free thiol chain (L_lig/2) in trans-NP is accounted for by a soft-repulsion term to account for the steric hindrance (labeled 3 in Fig. 2B). To model trans-NPs, we considered two different values of the surface-to-surface distance, d, between a pair of trans-NPs (Fig. S1, ESI†). At d = 0 (Fig. S1A, ESI†), the free thiol chain penetrates through the DDA/DDAB shell. Due to the steric hindrance from the shell, this conformation is assumed to be unstable. In the numerical calculation of the interparticle potential energy, this steric effect is incorporated as an excluded volume soft-repulsion for d = 0 to (L_lig/2). Beyond d > L_lig/2 (Fig. S1B, ESI†), this assumption is discarded.

Therefore, each of the cis/trans-NPs have an effective radius: R = R_NP + (L_lig/2), where R_NP is the radius of the NP core and L_lig is the functional length of the ligand. The arrangement of the azo-groups on the ligand-coated NP surface mimics the ‘points on a sphere’ model. Based on the calculations of dipole–dipole interactions and adhesion energy between two ligand-coated NPs, the number of ADT molecules on the NP surface is calculated to be approximately ∼19 ligands for the formation of reversible structures.⁶ We assume that the NPs are coated with 20 ligands, and the distribution of ADT on the NP surface is uniform.

2.2 Calculation of ligand-mediated effective interaction energy of a pair of ligand-coated NPs

The interactions between the ligands are non-specific in nature and a ligand on one nanoparticle can interact with multiple ligands on the other NP, if they are present within a favorable range of interaction. Therefore, the effective interaction energy for a pair of ligand-coated NPs is considered to be the sum of the interaction energies between each ligand–ligand pair on the mutually visible surfaces of the particles. Based on the placement of the ligands on cis/trans-NPs, the distance between a pair of ligands (r_ij) is calculated as the distance between the center of mass of the azo-groups of the ligands (Fig. 3). For an interacting pair of particles, mutually visible hemispheres of particles in a given particle-pair contribute to the interaction energy. For a single particle coated with 20 ligands, one hemisphere contains 10 ligands given the uniform distribution of ligands. Therefore, there are 100 interacting ligand–ligand pairs (N_pair) for a pair of particles, which contribute to the overall interaction potential. Each ligand–ligand pair interacts via the generalized LJ potential derived based on the DFT calculations.³⁶ We calculated the overall interaction potential for a pair of cis/trans-NPs using both analytical and discrete approaches.

Fig. 3 Schematic representations highlighting the distances between a pair of points (azo-groups) on the surfaces of a pair of (A) cis-NPs and, (B) trans-NPs. The schematic representations of the ligands show differences in orientations of the ligands in cis and trans conformations.

2.2.1 Discrete method. The discrete method considers the ‘points on a sphere’ model of the ligand-coated NP. Each point is considered to be the azo-group of the ligands. We modeled the placement of azo-groups on the ligand-coated NPs as a sphere with evenly distributed points, using the Fibonacci lattice method.⁴⁹Fig. 4A shows a pair of ligand-coated NPs depicted by spheres, with evenly distributed points on the surface. For each particle, the points on one hemisphere are colored green and on the other hemisphere as red. Each red point on particle 1 interacts with all green points on particle 2. The points are colored differently only to show the difference between the hemispheres. However, the nature of interactions is the same for all the points. The effective interaction energy [U(d)] for a pair of ligand-coated NPs 1 and 2 is therefore given by the sum of interactions between the ligands on the right hemisphere of particle 1, interacting with the ligands on the left hemisphere of particle 2. For N_lig/2 ligands on each interacting hemisphere, the total interaction can be written as:where,Here, [scr U, script letter U]_DFT(r_ij) is the interaction energy between a ligand i on particle 1 and a ligand j on particle 2, given by the generalized LJ potential from DFT calculations.³⁶ The range of [scr U, script letter U]_DFT(r_ij) determines whether a pair of ligands i and j on particles 1 and 2 would interact with each other. Also, r_ij is the distance between i and j, and d is the surface-to-surface distance between 1 and 2. A is a prefactor that determines the well-depth of the interaction, σ is an arbitrary length scale, and α and β are integers.³⁶

Fig. 4 (A) A schematic representation showing the discrete description of a pair of NPs. The red and green points represent the azo-groups on the right and left hemispheres of each particle respectively. Shown is the calculation of the distance between ligand pairs. The gray ring around the points represents the complete area of influence of the azo-groups on the NP surface. (B) A schematic representation of the volume within which the attached thiol chain of a ligand can move, making the effect of the azo-group spread over a wider area on the particle surface. (C) A schematic representation of the analytical description of a pair of NPs. The hemispherical surfaces S₁ and S₂ on particles 1 and 2, contributing to the effective interaction energy, are colored in darker shades. Shown is the distance between two arbitrarily chosen points P₁ and P₂ on S₁ and S₂, respectively, and their description in spherical coordinates.

Though the discrete method ideally maps the position of the azo-groups on the NP surface, it does not account for the movement of the ligand within a conical volume about the point of attachment on the surface, which allows the azo-groups to move over a fixed area over the surface of the NP. Fig. 4B shows a representation of the projection of the conical volume within which each ADT ligand is assumed to be free to move. Due to this effect, the azo-groups on the surface of the sphere can be considered to be displaceable within the vicinity of its original position, as shown by the gray ring around the points in Fig. 4A. Also, the ADT ligands are attached to the gold nanoparticle surface via thiol-bond linkages. The gold–thiol bond is a weak covalent bond due to which the ligands are free to move over the surface of the gold nanoparticle.^50,51 This information along with the flexibility of the ligand to move within the conical volume can be combined to hypothesize that the effect of the dipole moment of the azo-group can be considered to be continuous over the surface of the ligand-coated NP. In the discrete method of computing the interaction energy for a lower number of points (∼20), the positions of the points can affect the magnitude of the calculated interaction energy. Therefore, we considered the influence of the azo-group of the ligand to be continuous over the surface of the particle, in contrast with a point-like representation. This is implemented via the analytical approach for obtaining the interaction potential for a pair of NPs, as described further.

2.2.2 Analytical method. In the analytical method, we considered a uniform coverage of a specific number of ligands on the NPs (N_lig = 20). This uniform coverage is modeled as a uniform surface ligand density, ρ_lig. The overall interaction energy of a pair of cis/trans-NPs is then calculated by integrating over the interacting hemispheres in a ligand-coated NP pair. We considered the interactions between the exposed hemispheres S₁ and S₂, of a pair of ligand-coated NPs 1 and 2 (Fig. 4C). The derivation of the analytical expression describing the overall interaction energy between the particles is discussed further. The parameters used in the derivation are summarized in Table 1.

Table 1 A summary of parameters for deriving the analytical interaction energy between a pair of ligand-coated NPs

Parameter	Description
R NP	NP radius
L lig	Ligand/ADT length
	Effective radius for ligand-coated NP
d	Surface-to-surface distance between two ligand-coated NPs along the line connecting their centers
r = d + 2R	Center-to-center distance between two ligand-coated NPs
N lig	Number of ligands on the NP surface
	Ligand density on the NP surface
N pair	Number of interacting ligand pairs
r 12	Distance between points P1 and P2
	Dimensionless R
	Dimensionless d
	Dimensionless r
	Dimensionless r12

---

We first considered two points, P₁ on the surface S₁ of particle 1, and P₂ on the surface S₂ of particle 2. A spherical coordinate system is used (Fig. 4C).

P₁ = (R cos θ₁ sin ϕ₁,R sin θ₁ sin ϕ₁,R cos ϕ₁), θ₁ ∈ [0, 2π], ϕ₁ ∈ [0, π/2]

P₂ = (R cos θ₂ sin ϕ₂,R sin θ₂ sin ϕ₂,2R + d − R cos ϕ₂), θ₂ ∈ [0, 2π], ϕ₂ ∈ [0, π/2]

where, θ is the azimuthal angle and ϕ is the polar angle of a point on the ligand-coated NP. The distance (r₁₂) between points P₁ and P₂ is given by,The effective interaction energy between the ligand-coated NPs 1 and 2 is calculated as a function of their surface-to-surface distance (d) by integrating over the interacting hemispheres S₁ and S₂ according to the equation:In eqn (4), dS is given as:

dS = R² sin ϕdϕdθ (5)

For the hemispherical surface, θ ∈ [0,2π] and ϕ ∈ [0,π/2]. Therefore, the surface integral for the hemisphere is:

Thus, the effective interaction energy [U(d)] between the ligand-coated NPs can be written as,

Eqn (7) is solved numerically using the Mathematica software.⁵²U(d) is calculated for a range of d values until the magnitude of U(d) approaches zero. The equivalency of the discrete and analytical approach for [scr U, script letter U]_DFT(r₁₂) = 1 is demonstrated in ESI.†

2.3 Simulation setup

Coarse-grained Langevin dynamics simulations of cis/trans-NPs were performed in the NVT ensemble by using the HOOMD-blue software.⁵³ Each system consisted of N = 8000 particles. Each system was studied at different values of the reduced temperature, T* (ε_T/k_B units) (0.2, 0.15, 0.1, 0.08, 0.06, 0.05, 0.04, and 0.02) as well as for a range of values of the reduced number density, ρ* = N/L³ (σ_P⁻³ units) (1.25 × 10⁻⁴, 2.96 × 10⁻⁴, 1 × 10⁻³, and 8 × 10⁻³), modeled by varying the corresponding box size, L (σ_P units) (400, 300, 200, and 100). All simulations were performed by quenching the system from a high temperature, where the system is in a random configuration, to the specified temperature, T*. The number of NPs (N) is kept constant and the effect of density is studied by varying L, corresponding to a change in ρ*. Here, ε_T is the magnitude of the interaction energy, k_B is the Boltzmann constant, and σ_P = 1 is the diameter of the ligand-coated NP. We performed all simulations in cubic simulation domains with periodic boundary conditions applied in all directions. Each system was equilibrated for 2.5 × 10⁹ steps, followed by another 2.5 × 10⁹ production steps, during which the equilibrium properties of each system were monitored. A time-step of 0.0025 was used in all simulations.

2.4 Analyses techniques

The visual characterization of each system was done using the visual molecular dynamics (VMD) software.⁵⁴ The spatial arrangement and phase characteristics of each system were determined by calculating the radial distribution function (RDF) using eqn (8)The RDF gives the local density and the arrangement of particles with respect to the distance (r_ref) from a reference particle in comparison to the bulk density (ρ). The RDF is calculated by time-averaging over the last 20 000 frames recorded every 50 000 time-steps during the production cycle of each simulation. A bin size of 0.2σ_P was chosen for computing the distributions.

The average cluster size (M_n) of the colloidal aggregates, which refers to the typical or mean size of the individual clusters or groups formed by colloidal particles when they aggregate or coalesce, was also calculated. Particles are assigned to clusters if they are within a center-to-center distance of 1.5σ apart. This choice of cut-off was made because at 1.5σ, the potential is about the half of its minimum value. The average cluster size is then measured using:

where m_c is the number of clusters, i(k) is the number of particles in the kth cluster. The time averaging of M_n is done over the production cycle.

The compactness of the clusters is measured using the mean squared radius of gyration (R²_g), which provides a way for quantifying how the mass or components of an object are distributed relative to its center of mass. It is calculated as the mean squared distance between the constituent particles in a cluster.

where N_c is the number of particles in a cluster, i and j are two particles in the same cluster and r_i and r_i are the coordinates of i and j.

3 Results and discussion

3.1 Effect of functional ligand length

The cis and trans isomers of ADT vary in their spatial arrangements (Fig. 1A). Due to this difference, the free thiol chain is neglected in the theoretical description of cis-NP, but is introduced for the trans-NP as an excluded volume factor in the range d = 0 to L_lig/2 (orange region in Fig. 2B-3). The choice of the functional ligand length, which depends on the structural conformation of the ligand assuming that the thiol chains are flexible, also determines the effective radius of the ligand-coated NP. The length of ADT used in experiments is calculated to be 40 Å in the trans conformation, with thiol chains extended completely. However, since the alkyl chains are flexible they can undergo bending in the bulk solvent (Fig. S2A, ESI†). Therefore, we assumed that the ligand length can be reduced up to 10 Å, since the length of the azo-group in the linear trans-ADT molecule is ∼9 Å.^44,55 DFT calculations for a pair of ADT ligands showed an attractive interaction energy well for the trans-ADT.³⁶ However, the self-assembly studies of ligand-coated NPs in experiments do not show a significant assembly for the trans-NPs.⁶ Therefore, we hypothesized that the functional ligand length that gives a near zero attraction between a pair of trans-NPs is an ideal value for the model. We chose 25 Å to be the functional length of the ligand in the theoretical model (Fig. S2B, ESI†), and we demonstrated good qualitative agreement with experimental findings. Though a small attractive region is observed in this case at a distance of 10 Å, the magnitude of this attraction is very small (∼0.002k_BT) in comparison with the well-depth of the cis interaction (∼0.6k_BT) and can be neglected. Therefore, the optimal functional ligand length is chosen to be 25 Å in our model, which is between the calculated maximum length 40 Å and the assumed possible shortest length 10 Å of the ligand.

3.2 Effective potential energy for a pair of ligand-coated NPs

Using the analytical method (see Section 2.2.2), we calculated the effective interaction energy for a pair of cis/trans-NPs as a function of the surface-to-surface distance, d (Fig. 5A). The potential energy function obtained for the cis case (pink trace) resembles an LJ like form, where the particles do not interact at large distances (zero energy). As the particles come closer, they experience an attraction (negative energy) until a certain point, beyond which the excluded volume repulsion applies (positive energy). Using theoretical calculations, Klajn et al.⁶ reported that the minimum distance between two cis-ADT dipoles is ∼5 Å. In agreement with this, a minimum in the potential energy function is also obtained at ∼5 Å in our model. For the trans case (cyan trace), the potential energy function resembles the Weeks–Chandler–Andersen (WCA) form, which is a variation of the soft-sphere potential. As discussed previously, the model enforces an excluded volume repulsion due to the thiol chains for the trans-NP interaction in the range d = 0 to (L_lig/2) (12.5 Å). However, the repulsion is seen to take effect only around 6 Å (Fig. 5A; cyan trace). This is because for d > 6 Å, only a small number of ligand-pairs that are in the interaction vicinity (along the equator of the particle) contribute to the overall interaction. As the particles move closer (d < 6 Å), more ligand-pairs start contributing to the total interaction. This also makes the potential energy softer.

Fig. 5 (A) Effective interaction energy for a pair of cis/trans-NPs calculated using the analytical method (eqn (7)). (B) A comparison of the effective interaction energy for a pair of cis/trans-NPs calculated using the analytical method (bold lines), the potential energy expressions derived through a fit (labeled with +), as shown in eqn (11) and (12), and the expanded-Mie interaction potential in the HOOMD-blue software (dashed lines), given by eqn (14).

3.3 Model fit of the effective interaction energy

For the usage of the computed interaction energies in simulation studies, we obtained analytical expressions through a fitting procedure. The nonlinear least-squares Marquardt–Levenberg algorithm in the gnuplot software was used to perform the fit.⁵⁶ Based on the characteristics of the interaction energy traces, the cis-NP pair-energy was fit to a potential energy function of the LJ form (eqn (11)).

The trans-NP pair-energy was fit to an expression of the WCA form (eqn (12)). The trans-NP energy can also be fit to a soft-sphere like expression. We chose the WCA-like expression based on the pair potential-energy models available in the HOOMD-blue software.

where, C is a constant based on the values of n and m and is given as:The model parameters ε, σ, n, m for both cis and trans cases were obtained by fitting the calculated energy (from eqn (7); Fig. 5A) to the potential energy functions (eqn (11) and (12)) and are different from those obtained from DFT calculations (eqn (2)).

The model fit parameters for the cis/trans-NPs are given in Table 2. Due to the extended thiol-chains that are accounted for in the trans-NP model, the σ value for the trans-NP is higher than the cis-NP. This also accounts for the n and m values being smaller for the trans-NP. In the cis-NP, the excluded volume repulsion comes into effect at the surface of the ligand-coated NP, i.e., at the outer edge of the DDA/DDAB shell (gray region labeled 2 in Fig. 2A). Since the surface of the cis-NP is considered to be non-penetrable, eqn (11) goes to inifinity rather steeply. On the other hand, the excluded volume repulsion in the trans-NP comes into effect from a distance away from the surface of the ligand-coated NP (orange region labeled 3 in Fig. 2B). This orange region between the border of the thiol-chain shell and the DDAB shell is less crowded as it consists only of the extended thiol-chains of the ligands. Therefore, the steric hindrance of overlap is lesser, making eqn (12) softer, leading to smaller n and m values.

Table 2 The parameters for the interaction energy equations shown in eqn (11) and (12), for cis and trans NPs, respectively

	ε (kBT)	σ (Å)	σ*	n	m
cis	3.5	4.3	0.172	9	6
trans	0.05	6.75	0.27	6	4

---

3.4 Application of potential energy functions in the simulation model

The effective interaction energy calculated using eqn (7) and the parameters obtained by fitting this energy to eqn (11) and (12) are obtained in real units of energy (k_BT) and distance (Å). For the purpose of using eqn (11) and (12) in simulations, the model parameters (ε, σ) should be non-dimensionalized (ε*, σ*) with respect to the chosen units for energy and length. The unit of energy in simulations is chosen to be k_BT. Since ε is obtained in units of k_BT, the numerical value of ε* remains the same as ε in simulations. The unit of length in simulations is chosen to be the length of the ligand, L_lig = 25 Å. Therefore, the σ values are scaled down by a factor of 25 to σ* in simulations. The values of σ* are also reported in Table 2.

In HOOMD-blue, the expanded-Mie pair potential-energy expression (eqn (14)) is used to implement the eqn (11) and (12) in simulations.

The expanded-Mie potential shifts the potential to the surface of the particle using the parameter Δ. The value of Δ is set to the dimensionless diameter of the ligand-coated NP (2R* = 3.24) to shift the potential to the particle surface in our simulations. The value of r_cut is set to 4 for the cis case and 3.58 for the trans case where the energy approaches zero, for all interacting pairs. Fig. 5B shows a comparison highlighting the consistency between the effective potential energy as a function of the dimensionless surface-to-surface distance, d* (eqn (7)) (bold lines), the energy calculated using eqn (11) and (12) with the dimensionless fit parameters ε* and σ* (+ symbols), and the expanded-Mie potential-energy equation used in HOOMD-blue (eqn (14)) (dashed lines).

3.5 Self-assembled morphologies

We performed Langevin dynamics simulations for the cis/trans-NP systems separately for multiple values of T* and ρ* or L, as discussed in Section 2.3. The morphologies obtained from simulations were analyzed via the RDFs, the average cluster size, and the R_g of the clusters. Fig. 6 shows the particle configurations for a system of cis-NPs at selected values of T* and ρ*. The configurations for all conditions of T* and ρ* studied are shown in Fig. S3 (ESI†). In the cis-NP system, we observed a transition from a disordered phase at high T* to self-assembled aggregates at low T*, at all values of ρ*. At lower densities (ρ* = 1.25 × 10⁻⁴, 2.96 × 10⁻⁴), the aggregates are smaller. With increase in density (ρ* = 1 × 10⁻³), two or more aggregates diffuse and combine with each other to grow in size. The self-assembled morphologies obtained are dynamically arrested due to the quenching from a high temperature to a low temperature, leading to the formation of elongated clusters. At moderate T*, the kinetic energy is just sufficient to cause rearrangement of the particles within the clusters. At this condition, the elongated aggregates can rearrange to become more spherical with longer simulation time. In the trans-NP system, we observed a disordered phase at all conditions of T* and ρ*.

Fig. 6 Snapshots from simulations of the cis-NP system as a function of reduced temperature, T* and reduced density, ρ* (fixed N and varying L).

The RDF calculations of the cis-NP system are performed to understand the local arrangement of particles in each phase. Fig. 7A shows the RDF traces for the system at ρ* = 2.96 × 10⁻⁴ at multiple values of T*. The RDF for the disordered gas-like phase shows a single peak that quickly normalizes to bulk density (g(r) = 1). The system is disordered for T* = 0.15 (violet trace), suggesting that the thermal energy is sufficiently high to overcome particle interactions. The aggregated phases on the other hand showed multiple peaks before it normalizes to the bulk density. The change in temperature also affects the coarsening of the clusters which is also captured in the RDF traces. At T* = 0.1, the system is at the onset of aggregation. The clusters formed at this condition are larger in comparison to the clusters at lower T*, and have a well-ordered local arrangement (green trace) due to sufficient kinetic energy that permits local rearrangement of the particles. The RDF of this phase indicates well-ordered peaks spanning a longer distance before reaching the bulk density.

Fig. 7 (A) Shown is the plot of the radial distribution function [g(r)] of the system at ρ* = 2.96 × 10⁻⁴ at different values of T*. The simulation snapshots represent the system morphology for different traces in the plot. (1) Large ordered clusters, (2) small elongated clusters, and (3) disordered phase. Shown are the plots for (B) the average cluster size M_n, and (C) the mean squared radius of gyration R²_g of the system as a function of T* at different values of ρ* (fixed N, varying L).

The rapid convergence of the RDF trace to bulk density at T* = 0.06 and 0.08 (yellow and cyan traces) indicates that the clusters formed are smaller in size. These clusters are amorphous and elongated because the attraction between the particles are significantly stronger at lower temperatures and the kinetic energy is not sufficient to permit particle rearrangements. It is interesting to note that such amorphous precipitates are also observed in experiments at a higher methanol content, thereby making the inter-particle interactions stronger. Methanol being polar induces a higher dipole moment in cis-ADT which results in a stronger inter-particle attraction.⁶ The RDF traces for the other ρ* values follow the same trend (Fig. S4, ESI†).

Fig. 7B shows the average cluster size of the system as a function of T* at different values of ρ*. The average cluster size of the system increases with density at all values of T*. At low T*, the attraction between the particles is high, leading to the formation of self-assembled clusters. For lower densities (ρ* = 1.25 × 10⁻⁴, 2.96 × 10⁻⁴) at very low T*, the system forms clusters which are kinetically arrested and the growth is limited by slow diffusion. With increase in T*, the clusters gain more kinetic energy to diffuse and merge with other smaller clusters, causing an increase in the cluster size. The error bars show a moderate distribution of cluster size about the mean cluster size. At higher densities (ρ* = 1 × 10⁻³, 8 × 10⁻³), the system gets percolated forming a larger amorphous network. At high T* (0.2–0.15), the particles are disassembled and therefore show a smaller or zero cluster size. The error bars in this case show a large variation. The M_n traces therefore show a non-monotonic trend with T* for medium-to-low densities, which is more pronounced at ρ* = 2.96 × 10⁻⁴. The non-monotonic trend of the cluster size has been reported in previous studies.³⁷ We make sure that the cluster size has equilibrated for the systems at very low temperature by plotting the average cluster size with respect to time steps (Fig. S5, ESI†). We also study the distribution of the cluster size [P(n) vs. n] at lower densities, where P(n) is the probability of finding a cluster of size n (Fig. S6, ESI†). The distribution is done for a bin-size of 100 (clusters of size between 0–100, 100–200 etc. are grouped together). The distribution is right-skewed for lower densities, but does not provide further insights into the nature of the clusters.

A trend similar to M_n is observed for R²_g (Fig. 7C). A higher value of R²_g indicates a bigger or a more elongated cluster. At higher densities (ρ* = 8 × 10⁻³, 1 × 10⁻³), the system reaches a percolated state and therefore, the R²_g attains a constant value with a decrease in T*. However, with a decrease in density, R²_g starts showing a non-monotonic trend as a function of T*. The R²_g at ρ* = 1 × 10⁻³ and 2.96 × 10⁻⁴ are slightly higher than that at ρ* = 8 × 10⁻³. This is because the clusters are more elongated owing to a larger volume which the clusters can occupy. At a low density ρ* = 1.25 × 10⁻⁴, the clusters are smaller and also vary in their spatial arrangement due to the large volume available. The error bars of the R²_g trace (purple trace) showed higher values, indicating a broader distribution of R²_g. Overall, our model captured the important elements observed in the self-assembly studies of the ligand-coated NP system. The cis-NPs undergo aggregation while the trans-NPs remain unaggregated. The clusters formed by the cis-NPs follow a non-monotonic trend, as reported previously.³⁷ This model would be suitable to perform future studies related to the phase-behavior and also to validate the aggregation kinetics.

Since the quenching simulations lead to kinetically arrested morphologies, we performed additional simulations using the simulated annealing method, in order to understand the system behavior at equilibrium. In this procedure, the system is slowly cooled from a high temperature to the specified temperature. Snapshots of the system morphology obtained with the simulated annealing method are shown in Fig. S7 (ESI†). Compact quasi-spherical clusters were observed at a lower density (ρ* = 1.25 × 10⁻⁴), while the clusters diffuse and merge with other quasi spherical clusters at moderate densities (ρ* = 2.96 × 10⁻⁴, 1 × 10⁻³). At a higher density (ρ* = 8 × 10⁻³), a close packed crystal with long-range order is obtained. However, we note that quenching simulations are preferred over simulated annealing because they resemble a more realistic experimental set up which leads to the formation of amorphous clusters.

It is also important to point out that studies have been conducted on various other ligand-coated colloidal particle systems. One such example is a system of DNA-coated colloids, where the surface of the colloid is coated with mobile single-stranded DNA (ssDNA) ligands.⁵⁷ Such systems exhibit open or closed packed crystalline structures.^58,59 In this system, the ligands link via very specific hybridization between complementary ssDNA sequences. Therefore, if a ligand on one particle binds to another ligand on another particle, these ligands are then not available for further binding. The overall potential energy for the system therefore involves a combinational entropy term. On the other hand, in our system the ADT ligands interact via non-specific interactions, which is different from the ssDNA coated nanoparticles. Due to these nonspecific interactions, a ligand on one NP can interact with multiple ligands on the other NP, if they are present within a favorable range of interaction defined by the DFT potential (eqn (2)).

What is also critical in our work is the switching of the interaction when the ligands undergo trans-to-cis isomerization upon UV exposure. The main highlight of our work is to use quantum chemistry derived empirical potentials to model the change in morphology of the ADT-coated NPs during this transition (Fig. 2), which also manifests as a change in the interaction, thus providing predictive capabilities to our model. We conclude that the higher excluded volume of trans-ADT coated NPs when compared to cis-ADT coated NPs is responsible for the change in behavior during transition from visible to UV light.

4 Conclusions

LISA is a non-invasive method that can be applied to photosensitive-ligand-coated NPs to direct the reversible formation of morphologies of different size and spatial arrangements. However, to fine tune the parameters controlling LISA, it is crucial to understand the inter-particle interactions of ligand-coated NPs and to derive its mathematical form which can be applied in theoretical and computational studies. Here, we report the modeling of the ligand-mediated interaction energy between a pair of gold nanoparticles coated with photoresponsive ADT ligands. The model described several aspects of the nanoparticle–ligand system such as the arrangement of ligands on the NP surface, the effect of ligand conformation (cis/trans) and the functional ligand length on the total interaction energy. We followed discrete as well as analytical descriptions of the ligand arrangement on the NP and discussed its relevance to the model. The discrete model considered a point like arrangement of the azo-groups on the NP surface. To further incorporate ligand flexibility, an analytical description of the ligand arrangement on the NP was adopted. The calculated interaction energy was fit to generalized expressions for utility in conducting simulations. The interaction energy has an LJ-like form for the cis-NP, while it is a WCA-like form for the trans-NP. The model incorporated important experimental parameters like the range of interaction and excluded volume, which were otherwise not captured by the general potential expressions used in previous studies of LISA in ligand-coated NPs. The dimensionless fit expressions were applied in Langevin dynamics simulations to study the phase behavior of the system. The cis-NP system showed disordered as well as aggregated phases based on T* and ρ* values. The trans-NP system remained unaggregated at all conditions of T* and ρ*. The model captured similar trends with respect to self-assembled morphologies and the average cluster size, as observed in previous studies. Overall, our model captures the inter-particle interactions between ligand-coated NPs in a rigorous and realistic manner, and can be tuned to study LISA in other ligand–nanoparticle combinations. The model also lays the foundation to establish an understanding of the kinetics of aggregation and dissociation during the light-induced reversible clustering that can be tuned by controlling the periodicity of light exposure.³⁷

Author contributions

R. A. M. K.: conceptualization, data curation, formal analysis, investigation, methodology, visualization, and writing – original draft. P. K. J. and H. V.: conceptualization, supervision, funding acquisition, project administration, resources, and writing – review & editing.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We gratefully acknowledge the financial and computational support provided by the National Science Foundation (NSF) EPSCoR award (OIA-1757371; H. V.). We also thank the support through Grant No. CRG/2021/001260 (P. K. J.) by the Science and Engineering Research Board, India.

Notes and references

1. M. A. C. Stuart, W. T. Huck, J. Genzer, M. Müller, C. Ober, M. Stamm, G. B. Sukhorukov, I. Szleifer, V. V. Tsukruk and M. Urban, et al. , Nat. Mater., 2010, 9, 101–113.
2. P. Theato, B. S. Sumerlin, R. K. O'Reilly and T. H. Epps III, Chem. Soc. Rev., 2013, 42, 7055–7056.
3. M. Mrinalini and S. Prasanthkumar, ChemPlusChem, 2019, 84, 1103–1121.
4. S. Municoy, M. I. Álvarez Echazú, P. E. Antezana, J. M. Galdopórpora, C. Olivetti, A. M. Mebert, M. L. Foglia, M. V. Tuttolomondo, G. S. Alvarez and J. G. Hardy, et al. , Int. J. Mol. Sci., 2020, 21, 4724.
5. Z. Li, Q. Fan and Y. Yin, Chem. Rev., 2021, 122, 4976–5067.
6. R. Klajn, K. J. Bishop and B. A. Grzybowski, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 10305–10309.
7. Y. Yan, J. I. Chen and D. S. Ginger, Nano Lett., 2012, 12, 2530–2536.
8. P. K. Kundu, D. Samanta, R. Leizrowice, B. Margulis, H. Zhao, M. Börner, T. Udayabhaskararao, D. Manna and R. Klajn, Nat. Chem., 2015, 7, 646–652.
9. D. Manna, T. Udayabhaskararao, H. Zhao and R. Klajn, Angew. Chem., Int. Ed., 2015, 127, 12571–12574.
10. Y. Liu, X. Han, L. He and Y. Yin, Angew. Chem., Int. Ed., 2012, 51, 6373–6377.
11. Nonappa, J. S. Haataja, J. V. Timonen, S. Malola, P. Engelhardt, N. Houbenov, M. Lahtinen, H. Häkkinen and O. Ikkala, Angew. Chem., Int. Ed., 2017, 129, 6573–6577.
12. Q. Zhang, W.-Z. Wang, J.-J. Yu, D.-H. Qu and H. Tian, Adv. Mater., 2016, 29, 1604948.
13. Y. Yu, D. Yu and C. A. Orme, Nano Lett., 2017, 17, 3862–3869.
14. M. Wang, L. He, W. Xu, X. Wang and Y. Yin, Angew. Chem., Int. Ed., 2015, 54, 7077–7081.
15. T. Bian, L. Shang, H. Yu, M. T. Perez, L.-Z. Wu, C.-H. Tung, Z. Nie, Z. Tang and T. Zhang, Adv. Mater., 2014, 26, 5613–5618.
16. D. Samanta and R. Klajn, Adv. Opt. Mater., 2016, 4, 1373–1377.
17. X. Yao, T. Li, J. Wang, X. Ma and H. Tian, Adv. Opt. Mater., 2016, 4, 1322–1349.
18. T. Bian, Z. Chu and R. Klajn, Adv. Mater., 2020, 32, 1905866.
19. M.-M. Russew and S. Hecht, Adv. Mater., 2010, 22, 3348–3360.
20. C. Garcia-Iriepa, M. Marazzi, L. M. Frutos and D. Sampedro, RSC Adv., 2013, 3, 6241–6266.
21. R. Klajn, P. J. Wesson, K. J. Bishop and B. A. Grzybowski, Angew. Chem., Int. Ed., 2009, 48, 7035–7039.
22. H.-B. Cheng, S. Zhang, J. Qi, X.-J. Liang and J. Yoon, Adv. Mater., 2021, 33, 2007290.
23. M. Di Martino, L. Sessa, R. Diana, S. Piotto and S. Concilio, Molecules, 2023, 28, 3712.
24. X. Chen, L. Lin, Z. Li and H.-B. Sun, Adv. Funct. Mater., 2022, 32, 2104649.
25. R. Klajn, K. J. Bishop, M. Fialkowski, M. Paszewski, C. J. Campbell, T. P. Gray and B. A. Grzybowski, Science, 2007, 316, 261–264.
26. S. Das, P. Ranjan, P. S. Maiti, G. Singh, G. Leitus and R. Klajn, Adv. Mater., 2013, 25, 422–426.
27. A. Zep, M. M. Wojcik, W. Lewandowski, K. Sitkowska, A. Prominski, J. Mieczkowski, D. Pociecha and E. Gorecka, Angew. Chem., Int. Ed., 2014, 53, 13725–13728.
28. F. A. Jerca, V. V. Jerca and R. Hoogenboom, Nat. Rev. Chem., 2022, 6, 51–69.
29. J. Li and X. Jia, Nanomaterials, 2023, 13, 2562.
30. M. Younis, S. Ahmad, A. Atiq, M. Amjad Farooq, M.-H. Huang and M. Abbas, Chem. Rec., 2023, 23, e202300126.
31. J. Karcher and Z. L. Pianowski, Chem. – Eur. J., 2018, 24, 11605–11610.
32. A. Köhntopp, A. Dabrowski, M. Malicki and F. Temps, Chem. Commun., 2014, 50, 10105–10107.
33. J.-W. Lee and R. Klajn, Chem. Commun., 2015, 51, 2036–2039.
34. B. G. van Ravensteijn, I. K. Voets, W. K. Kegel and R. Eelkema, Langmuir, 2020, 36, 10639–10656.
35. G. S. Hartley, Nature, 1937, 140, 281.
36. R. Singh and P. K. Jha, Comput. Theor. Chem., 2021, 1206, 113492.
37. P. K. Jha, V. Kuzovkov, B. A. Grzybowski and M. O. De La Cruz, Soft Matter, 2012, 8, 227–234.
38. R. Zhang, D. A. Walker, B. A. Grzybowski and M. Olvera de la Cruz, Angew. Chem., Int. Ed., 2014, 126, 177–181.
39. D. Pagan and J. Gunton, J. Chem. Phys., 2005, 122, 184515.
40. J. M. Ilnytskyi, A. Slyusarchuk and M. Saphiannikova, Macromolecules, 2016, 49, 9272–9282.
41. N. C. Craven, J. B. Gilmer, C. J. Spindel, A. Z. Summers, C. R. Iacovella and C. McCabe, J. Chem. Phys., 2021, 154, 034903.
42. P. Khan, R. Kaushik and A. Jayaraj, ACS Omega, 2022, 7, 47567–47586.
43. J. Henzl, M. Mehlhorn, H. Gawronski, K.-H. Rieder and K. Morgenstern, Angew. Chem., Int. Ed., 2006, 45, 603–606.
44. E. Merino and M. Ribagorda, Beilstein J. Org. Chem., 2012, 8, 1071–1090.
45. I. Conti, M. Garavelli and G. Orlandi, J. Am. Chem. Soc., 2008, 130, 5216–5230.
46. Y. Ootani, K. Satoh, A. Nakayama, T. Noro and T. Taketsugu, J. Chem. Phys., 2009, 131, 194306.
47. I. C. Merritt, D. Jacquemin and M. Vacher, Phys. Chem. Chem. Phys., 2021, 23, 19155–19165.
48. J. Fu, X. Zhang, B. Miao and D. Yan, J. Chem. Phys., 2017, 146, 164901.
49. Á. González, Math. Geosci., 2010, 42, 49–64.
50. H. Häkkinen, Nat. Chem., 2012, 4, 443–455.
51. Y. Xue, X. Li, H. Li and W. Zhang, Nat. Commun., 2014, 5, 1–9.
52. S. Wolfram, Mathematica: a system for doing mathematics by computer, Addison Wesley Longman Publishing Co., Inc., 1991.
53. J. A. Anderson, J. Glaser and S. C. Glotzer, Comput. Mater. Sci., 2020, 173, 109363.
54. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
55. H. Koshima, N. Ojima and H. Uchimoto, J. Am. Chem. Soc., 2009, 131, 6890–6891.
56. T. Williams, C. Kelley and many others, Gnuplot 4.4: an interactive plotting program, https://gnuplot.sourceforge.net/, 2010.
57. S. Angioletti-Uberti, P. Varilly, B. M. Mognetti and D. Frenkel, Phys. Rev. Lett., 2014, 113, 128303.
58. E. Auyeung, T. I. Li, A. J. Senesi, A. L. Schmucker, B. C. Pals, M. O. de La Cruz and C. A. Mirkin, Nature, 2014, 505, 73–77.
59. R. V. Thaner, Y. Kim, T. I. Li, R. J. Macfarlane, S. T. Nguyen, M. Olvera de la Cruz and C. A. Mirkin, Nano Lett., 2015, 15, 5545–5551.

---

Footnote

† Electronic supplementary information (ESI) available: Including supplementary methods used in the model, additional figures for clarity and additional snapshots from simulations of the cis and trans NP systems. See DOI: https://doi.org/10.1039/d4sm00459k

---