key: cord-0979167-24ffohyd
authors: Lieu, Uyen Tu; Yoshinaga, Natsuhiko
title: Inverse design of two-dimensional structure by self-assembly of patchy particles
date: 2021-09-21
journal: The Journal of chemical physics
DOI: 10.1063/5.0072234
sha: ac716ad5b1e1144dd810d44a67b7844e559af311
doc_id: 979167
cord_uid: 24ffohyd

We propose an optimisation method for the inverse structural design of self-assembly of anisotropic patchy particles. The anisotropic interaction can be expressed by the spherical harmonics of the surface pattern on a patchy particle, and thus arbitrary symmetry of the patch can be treated. The pairwise interaction potential includes several to-be-optimised parameters, which are the coefficient of each term in the spherical harmonics. We use the optimisation method based on the relative entropy approach and generate structures by Brownian Dynamics simulations. Our method successfully estimates the parameters in the potential for the target structures, such as square lattice, kagome lattice, and dodecagonal quasicrystal.

The self-assembly of nano-and colloidal particles is a spontaneous organisation of the small particles into structures [1] . The self-assembly forming complex patterns is of great interest because of the promising applications in materials engineering, such as photonics, energy storage devices, tunablerheology fluids [1, 2] . The special properties of these materials originate from the material structures and the physicalchemical properties of the components made of the materials [3, 4] . Understand, control and predict self-assembly structure from given building blocks is a challenging goal in the field of soft matter. The self-assembly of complex structures, particularly open structures such as kagome lattice [5] , honeycomb lattice [6] , diamond structure [7] , icosahedral quasicrystal [8] , however, has not fully understood yet. The main difficult issue is to determine building blocks which are capable of forming the desired structure, and the underlying kinetics and mechanism of the assemble process, i.e. the interaction of such building blocks and the conditions, e.g. temperature, dispersity will determine the outcome structure. The interaction of building blocks depends on the block and the solvent. The building block can be isotropic [8, 9] or anisotropic [10] . In the later case, the building block is varied in shape [11] , local sphere clusters [12] , patchy particle [1, 13] , etc. The patchy particle is often described as a spherical particle patterned with anisotropic surface, or attached with interacting patch on its surface [1, 13, 14] . Recent developments in synthesis and fabrication techniques have enabled the realisation of those patchy particles [5, 15] . The assembly of patchy particles has led to the discovery of new order structures, phase diagram, and the extraction of some general features such as the formation of hierarchical assembly [16] . Discovering the relation of particle design and selfassembled structure employs forward approach or inverse approach. In the conventional 'forward' self-assembly approach, a specific potential or anisotropic particle is used to discover the assembled structure. Conversely, in 'inverse' approach a class of computational technique is used to discover the suitable type of potential or anisotropic particle which is capable of forming a desired target structure [4, 17, 18] . For patchy particles, the infinite degrees of freedom of patchiness result in a great flexibility to achieve several complex structures. However, it also poses questions on the inverse problem: how to choose a suitable particle design for a specific structure among several possible particle designs. To determine suitable interactions for stabilising self-organisation, the optimasation techniques have been continuously developed for both isotropic interaction [6, [17] [18] [19] [20] [21] [22] [23] and anisotropic interaction [24] [25] [26] [27] .

In this study, we have proposed an inverse design strategy for the self-assembly of monodispersed patchy particles by Brownian dynamics simulation. The algorithm is based on relative entropy method [17] to optimise the yield of the target structure in the design parameter space. We apply the model to find suitable patchy particle designs for given target structures. The patchy particle in the study is a spherical particle whose surface pattern is based on spherical harmonics, thus any pattern can be described by a linear combination. [14] Such pattern consists of two type of regions which is an analogical to a polarised surface or even the surface charges of complex biological molecules [13, 28, 29] . Then one can set either repulsive or attractive interaction to the like/opposite patches, for example setting like patches repulsive and opposite patches attractive, which is similar to the charged interaction. The main feature of our patterned patchy particle system involves the complex interplay between the attractive and repulsive anisotropic interaction, and the capability of systematically exploiting the relation between the patchy particle symmetry and the self-assembly. Compared to the patchy particles with specific narrow attractive sites (sticky-patch model) where the patch number, patch arrangement and patch type are important features determining the assembled structure and properties, the particle in our approach are less involved with those factor. We expect that the experiment can be carried out more flexible by tuning the external parameters (solvent pH, salt concentration) [13, 15] and the fabrication required less strict conditions on such as a sophisticate geometry and chemical selectivity.

In this study, the pattern of patchy particle i is described by

(i) } includes the information of the local orientations of particle i [], and is the l-fold contraction of two l th rank tensors. Then for a pair of Y lm particles put at the distance r, the anisotropic interaction Ξ lm is described as

For simplicity, we assume that the potential u lm for a pair of particles Y lm can be decomposed into a distance-dependent term and an orientation-dependent term as u lm = u M (r)Ξ lm (r, Ω). Here u M (r) is similar to the Morse potential, and Ξ lm (r, Ω) is dependent on the mutual orientation Ω of the pair particle. We

In the inverse problem, the to-be-optimised potential includes the interactions of the several types of patchy particles:

where the parameter θ lm is the weight of each potential type to be optimised. For demonstration, the potential includes five candidates corresponding to five types of patchy particles Y 10 , Y 20 , Y 22 , Y 44 , Y 55 , as shown in Fig. 1 . The particle Y 10 and Y 20 can be referred as Janus particle and triblock particle, respectively. They are axisymectric about the polar axis (n (3) ). The particle Y 22 , Y 44 and Y 55 has twofold, fourfold, and fivefold symmetry in the equator plane spanned by (n (1) , n (2) ). The patchiness of a particle can be determined by its positive (red) and negative (blue) pattern. In this model, θ > 0 means that different colour patch is attractive and similar one is repulsive, while θ < 0 means similar colour patch attractive and different is repulsive. The sign of θ implies different interaction and physics of the self-assembly, therefore it can be considered as a prior knowledge for the design of patchy particle. Since the sign of θ is often fixed by materials of the patches, we can set either θ ≤ 0 or θ ≥ 0. We may also set the sign unconstrained. The detail of the anisotropic function Ξ lm , the Morse potential u M , and the Week-Chandler-Anderson potential u WCA preventing the overlapping of particle and can be found in Appendix. The self-assembly of patchy particles is performed by Brownian dynamics [30] under annealing. In each simulation, the spherical particles are confined to a flat plane while rotate freely in three dimensions. The periodic boundary condition is applied on the two-dimensional space. The number of particles in each iteration is N = 128 for square lattice and kagome lattice, N = 240 for dodecagonal quasicrystal. Initially the positions and orientations of particles are randomly distributed. Temperature is decresed from T max to T min with an interval of ∆T . The values of T max , T min , ∆T are respectively 0.5, 0.05, 0.01 for square lattice and kagome lattice, and 1, 0.2, 0.0125 for dodecagonal quasicrystal.

n (2) n (1) FIG. 1: The patchy particles described by spherical harmonic Y lm (top row) and the corresponding patchiness (middle and bottom rows). The orientation of a particle is characterised by the local orthonormal bases (n (1) , n (2) , n (3) ). The view views along the polar axis n (1) of each particle are shown in the bottom row.

Our optimisation process is based on the relative entropy approach [17] . According to Jadrich et al. [17] , a model for isotropic interaction is built as follows. A particle configuration R is described by a set of N particles' positions R = [r 1 , r 2 , ..., r N ]. Let P tgt (R) denote the probability distribution of a target structure, and P(R|θ) denote the probability distribution for realising the configuration R given θ. Here P(R|θ) follows Boltzmann distribution P(R|θ) = e −βU(R|θ) Z(θ) , in which the configurational partition function Z(θ) is a normalisation factor as Z(θ) = e −βU(R|θ) dR, and U(R|θ) is the tunable potential. One way to measure the distance between those two distributions is to use Kullback-Leibler divergence [31] . With the use of f (R) P(R|θ) = f (R)P(R|θ) dR, the Kullback-Leibler divergence from P(R|θ) to P tgt (R) can be written as D KL (P tgt ||P) = P tgt (R) ln P tgt (R) P(R|θ) dR = ln P tgt (R) P tgt (R) − ln P(R|θ) P tgt (R) .

(2)

In the optimisation scheme, θ is tuned to minimise D KL (P tgt ||P). The search for the local minima can be conducted by the gradient descent method. Then the next point of θ in the iteration is chosen by following the steepest descent so that:

where α is the parameter controlling how far the point moves along the gradient descent curve. In this study we empirically set α = 0.05. Applying of Boltzmann distribution and partition function for the gradient term in Eq. (3), we can rewrite the term as the following form

In this study, the interaction of the particle is dependent on both the position and orientation of the particle, therefore, information on the orientation of the particle is required for the optimisation. When orientation of the target structure is available, Eq. (4) is simply replaced by

As we discuss in Sec. II C, we also consider the target structure whose orientation is not measured. In this case, we treat the orientation of the target structure as a hidden variable. The optimisation is interpreted as minimisation between P tgt (R) and the probability distribution that is marginalised over the orientation, P(R|θ ) = P(R, n|θ )P(n). Then, we use, instead of Eq. (5),

where P(n) is estimated from the generated structures for each iteration.

The weight θ affects the energy scale and how the patches on the particle interact. For each data set of optimisation, θ is limited to the range of |θ | = [0, 1.2]. We also investigate the behaviour of θ when constraint on the sign θ ≥ 0 or θ ≤ 0 is applied.

In order to thoroughly evaluate the inverse design optimisation, the target structure is categorised into three groups in decreasing order of target's information: ground truth target, synthesised target, and minimal synthesised target. (i) The ground truth target is prepared by stabilising the particles using the potential in Eq. (1) with given values of θ. Although the annealing scheme is chosen slow enough, the assembled structure still has defects and thermal fluctuations. The target includes all information about the positions and orientations of particles. It should be noticed that the calculation of energy of the assemblies requires the information of the position and the orientation (i.e. three local orthonormal basesn (1) ,n (2) , n (3) , as shown in Fig. 1 ) of the particles. (ii) The synthesised target is prepared based on the ground truth, however, the particle position is perfectly set according to the unit cell of a targeted lattice, the particle surface distance is set zero; then the orientations of the particles are set corresponding to an energy minimising structure. (iii) In the final case, only the position of the synthesised target is given to further check the capability of our proposed model.

In detail, we choose square lattice, kagome lattice and a dodecagonal quasicrystal as target structures (Fig. 2) . The square lattice is a basic structure. The kagome lattice is an open structure which is difficult to create via self-assembly and it is supposed to have novel optical and mechanical properties [5] ; dodecagonal quasicrystal is also difficult in terms of their aperiodic nature [32] . The ground truths of these structures, namely the square lattice, kagome lattice, dodecagonal quasicrystal can be obtained by setting the non-zero θ in Eq.

(1) as θ 10 = 1, θ 20 = −1, θ 55 = −1, respectively. As for the synsthesised targets, particle positions and orientations of the square and the kagome lattice are perfectly assigned based on their unit cells. The synthesised dodecagonal quasicrystal target is constructed as an approximant crystal by suitable packed dodecagonal motif [33] , however the orientation is too difficult to perfectly assign. We overcome this problem by numerically stabilising the approximant crystal under the potential type same to the one used for ground truth. Finally, the minimal targets, which have only information of particle position, are also considered. In this case, the orientation of target is a hidden parameter. However, because the orientation of the target is required for each iteration, we statistically interpolate the orientation fields from the generated structure and assign to the target. The algorithm of the interpolation of orientation can be found in Appendix D. FIG. 2: Illustration of target structures include square lattice, kagome lattice, dodecagonal quasicrystal considering ground truth (left), synthesised structure (middle), and minimal structrure (right). A Fourier transform of the ground truth dodecagonal quasicrystal with 12-fold symmetry is given inset. Other insets are enlarged structures. Particle colours are coded with number of nearest neighbours, the blue and red arrows are the orientations of the particles in the local direction n (3) and n (1) (Fig. 1) . The voronoi tesselation is also embeded.

The optimisation starts with a random set of the to-beoptimised weight θ, followed by a dynamical self-assemble process. Then the structural data at the final timestep, which is called generated structure, will be used for Eq. (3) to update the θ for the next iteration. The number of iterations is set around 20, 40 or hundreds. To compare the generated structures and the target, we generate the structure with the estimated parameters and then check the local structure, the number of neighbours and the radial distribution function.

The performance of the optimisation method is confirmed by inverse design of the ground truth target structures. Such ground truth target structures are prepared from a specific Y lm potential as given in Eq. (1) . The information of the ground truth target includes both positions and orientations of all particles. For the ground truth square lattice target assembled from Y 10 patterned particle, Fig. 3 shows the typical behaviours of the constrained θ ≥ 0. As shown in Fig.  3(a,c) , a solution containing the only non-zero θ 10 ≈ 1 is quickly obtained, and the generated structure is identical to the ground truth target. We also prepare a synthesised lattice whose surface distance between particle is set zero, and orientation of polar basis n (3) is exact. Compared to the ground truth case, the behaviour of θ in Fig. 3 (d) is qualitatively similar. However the weight is θ 10 = 1.2, which is the upper limit of θ. The reason of this difference can be explained as follows: according to the estimation of θ given in Eq. (3) and (4), the terms ∇ θ U of the generated structure and the target is distance-dependent. The distance between particles in the synthesised target is smaller than that of the ground truth (the first peak of radial distribution function for synthesised target and ground truth is 2a and 2.08a respectively; data are not shown). Since the interaction decays with distance, the term ∇ θ U tgt of synthesised target has greater magnitude as shown in Fig. 3(b,e) , leading to ∇ θ 10 U − ∇ θ 10 U tgt = 0 for the synthesised lattice optimisation. As a result, θ 10 keep increasing until it reaches the upper limit. The result for the synthesised target suggests that the optimisation scheme can work when the distance between particle in the target is not strictly assigned, i.e. similar behaviour of θ is expected when the optimisation is conducted for a synthesised target whose distance between particle is larger or smaller. Additionally, it should be noted that in synthesised target, we assign only the polar basis (the n (3) in Fig.1 ) of particle while the information of azimuthal bases (n (1) and n (2) direction) are random; because the Y 10 particle is axisymetric, so the information of the polar basis is enough to attain a minimum energy structure. This result suggests that the model can work when the target contains partial information of the orientation.

The result for kagome lattice is given in Fig. 4 when the sign of θ is arbitrarily assigned. The ground truth target is obtained by θ 20 = −1. The optimisation is capable of finding the correct solution whose nonzero value θ 20 ≈ −0.9 for ground truth and θ 20 = −1.2 (lower limit) for synthesised target. The difference between the value of θ 20 in such two cases is mostly due to the particle distance in the two targets, as mentioned in previous paragraph. One can see that the generated structure contains quite a lot of defects, perhaps this is the reason why |θ 20 | for ground truth is slightly smaller than the expected value.

The estimation for dodecagonal quasicrystal is given in Fig.  5 . The ground truth is obtained by considering the potential with θ 55 = −1. For synthesised target, the same potential is used to numerically stabilise a position-fixed approximant crystal. Different from the synthesied target of square lattice and kagome lattice, we do this process to find the minimum energy state of the approximant because it is too difficult to assign the orientations to the particles. When θ ≤ 0, after a few iterations the optimisation model is capable of finding the suitable potential with θ 55 ≈ −0.9 for ground truth and θ 55 ≈ −0.7 for synthesised target. Although these values deviate from the expected value around −1, the convergence of θ behaves similarly and the only nonzero parameter for the estimation refers to the Y 55 patchy particle design. It is confirmed that the generated structures by those parameters have features of a dodecagonal quasicrystal, which are the 12-fold symmetry in their Fourier transforms in Fig. 5(b,d) , and the similarity with the radial distribution of the ground truth target as in Fig. 5(e) . Regarding the difference in estimation result, one possible reason is that the number of particles used in those simulation is quite small, therefore the periodic boundary condition can inhibit the growth of the quasicrystal and the relaxation of defects. For the synthesised target, the difference in the obtained value and expected value is more obvious. As mentioned previously, such effect is also caused by the particle distance in the generated structure and the target, which is also illustrated in Fig. 5(e) .

The defects in a generated structure have certain effect on the estimation of θ. The defect is manifested in the form of grain boundary and point defect. For the square lattice case, the generated structure at the latter iterations contains a few defects and the obtained solution θ 10 ≈ 1 (θ ≥ 0) is almost similar to the potential exclusively used to create the ground truth target. For the kagome lattice, it can be seen that the generated structure have quite a lot of defects. The generated structure of quasicrystal is not as clear as the ground truth, possibly because of the small number of particles. In these two case, the obtained solution (θ 20 for kagome lattice and θ 55 for dodecagonal quasicrystal) have smaller |θ | than expected. There is a trade-off between accuracy and computational cost, one may adjust the annealing setting of the simulation to reduce defect, however, the simulation cost increases. In this study, although the number of particles seems small and the annealing is quite fast, the results show sufficient ingredients for designing the particle.

For the minimal target, the knowledge on possible orientation of particle is unknown. Therefore we extract the orientation field of the generated structure and pass this information to the target. Figure 6 illustrates the behaviours of θ and the possible solutions for square lattice and kagome lattice target. For the square lattice, when the constraint θ ≥ 0 is applied (different patch attractive condition), we obtain two different sets of pa- (c,f) The snapshots are generated using the parameters at the initial condition and the estimated result. The colour of particle is the number of nearest neighbours of each particle. The blue arrows are the orientations of the particles in the local polar direction n (3) (Fig. 1 ).

rameters. The nonzero parameters in these sets are θ 10 = 1.2 and θ 22 = 1.2, respectively. It suggests that aside from the Y 10 type as estimated in Fig. 3(a) , the other type of patchy particle, Y 22 , can assemble into a square lattice. This patchy particle type composes of four patches around equator, and the patches with different color are attractive (see Fig. 1 ). The capability of forming a square lattice for such kind of particle is intuitively understandable. In the case of kagome lattice target with the constraint θ < 0, as depicted Fig. 6(c) , we obtain the unique estimation in which there are two nonzero parameters as θ 10 = θ 20 = −1.2. The result suggests that the particle design is a combination of particle type Y 10 and Y 20 . This result is different from the solution of Y 20 found in previous sections (Fig. 4) .

In order to evaluate the contribution of each type, we have independently conducted three groups of simulation using the particle patterned with Y 10 , Y 20 , and Y 10 +Y 20 . Figure 7 shows a quantitative comparison of the structures obtained by the three particle types, by analysing the distribution of the number of nearest neighbours and the six-fold bond-orientational order parameter [41] |ψ 6 |. If a particle is a part of kagome lattice, it has four nearest neighbours and |ψ 6 | = 1. From the figure, for three types of patchy particle, the significant rate of n nb = 4 and |ψ 6 | ≥ 0.7 suggests that the structure is a kagome lattice [5] . In the case Y 10 + Y 20 , compared with the particles patterned only with Y 10 or with Y 20 , the number of more-than-4-neighbour particles smaller while the number of less-than-4-neighbour particles larger, implying that the structure has more open local structure. Moreover, the system has more particles whose ψ 6 > 0.9 , suggesting the kagome lattice is clearer. This result reveals that the combination of Y 10 and Y 20 is able to create a much better kagome lattice than using only one of them. As depicted in Fig. 7 , under the combined potential, the resulted kagome lattice has less defects, more open structure, and the shape of the unit cell is clearer than the case using the single component Y 10 or Y 20 . A qualitatively similar phenomenon is observed for the enhanced self-assembly kagome lattice from tri-block particle whose the self-propel activity is added to the pole direction of tri-block particle [34] . The tri-block particle and the adding of self-propulsion can be respectively referred to the Y 20 and Y 10 in our study. Although the nature of enhancement by the self-propulsion [34] is different from the patchy particle Y 10 , they are similar in terms of symmetry of the interaction between the particles. These results suggest that our study is capable of finding new designs for patchy particle so that the self-assembly is enhanced.

In summary, we have developed a relative entropy-based method for inverse structural design of patchy particle for a given target structure. The type of a patchy particle is described by the spherical harmonic. The pairwise interaction potential includes several to-be-optimised parameters, which are the coefficient of each term in the spherical harmonics. We successfully estimate patches necessary to reproduce the targets such as two-dimensional square lattice, kagome lattice and dodecagonal quasicrystal. The method also works for hidden information for the square and kagome lattices, i.e. the target contains only the minimal information of the position of particles. The estimation is dependent on the type of a target structure, namely, prior knowledge of the structure that we want to reproduce. The prior knowledge of the range of the parameters is also crucial for the estimation. Here, we discuss those effects, as well as the possible extension of the model.

We have considered three kinds of target structure: the ground truth target comprises both the particle position and orientation, the synthesised target has its particle position set in perfect order and orientation obtained via minimum energy, and the minimal target has only information of the particle position. For the periodic target like square lattice and kagome lattice, the optimisation model is capable of finding the suitable solution design so that the desired targets can be obtained. Regarding the performance of the estimation for three types of target, it shows that as the more information on the target is provided, the faster the optimisation process is. For example in the case of square lattice target (Fig. 3 and Fig. 6 ), for the ground truth and synthesised targets, the stable and converging solution can be consistently obtained after just a few iterations, whereas the minimal target requires more iterations. The minimum target does not have data on the orientation, so the optimisation involves searching on a higher dimension space. As a result, it is more difficult to estimate the parameters. We also observe that the occurrence for finding the correct solution of minimum target is often lower than the other target. Among many independent estimation processes for each square lattice target, the solutions for ground truth and synthesised ones are all consistent and can reproduce a square lattice structure, while for the minimum target cases about 50% of the estimations can. However, such disadvantage is compensated by the fact that more possible solutions and "better" solutions can be found. In detail, when the mini- mum targets are considered, we have found two possible particle designs for the square lattice, three for the kagome lattice in which a combination design improves the lattice. On that account, one can consider different levels of prior knowledge for target structure. Here is an example for the self-assembly of kagome lattice. From the lattice, one can think that there can be at least two ways to design the particle's patches and the corresponding orientation of those particle in the lattice. Intuitively, the design based on narrow sticky patch consists of four small patches located on the same plane at alternating interval of 60 • and 120 • so that they fit the four contact points [35] [36] [37] [38] . However, self-assembly of those particles shows that a rhombus lattice is more favourable and formed instead of kagome one [35] [36] [37] [38] . The kagome lattice structure is formed only when constraint on patch selectivity is added [37, 38] . Another solution is the axisymmetric particle, then the orientation of particle can be inferred, which is similar to the one in synthesised target. Using such kind of information for the target will result in the finding of Y 20 particle design, in which two opposite identical caps are attractive while the equatorial part is repulsive [5, 34, 39] .

In the case of minimum target dodecagonal quasicrystal, the optimisition is not able to find the expected solution. Aside from the difficulty caused by the hidden orientation as discussed above, there are several factors hindering the process. The target does not have information of particle orientation, so it is necessary to assign the orientation from the generated structure to the target. In the generated structure, the effect of thermal fluctuation is inevitable. Such fluctuation is also captured and passed to the target. In other words, simply imagine that the number of patches on the particle Y 55 of quasicrystal is large, the size of the patch is thus small, and then the fluctuation may lead to a considerable effect. So the effect of this fluctuation on the estimation outcome is more severe than that of the square lattice and kagome lattice. Another factor is that compared to the formation of the aperiodic lattice like dodecagonal quasicrystal is difficult and thus requires more investigation on the suitable temperature range, annealing rate, as well as the influence of the other potential components.

In this study, the sign of θ involves the way the patch interacts, thus the potential type and fabrication technique. The patchy particle has two kind of regions (red and blue) on its surface, and different-kind region is set attractive for θ > 0, while same-kind region is set attractive for θ < 0. It should be noted that the self-assembly is completely different when the sign changes. Applying constraint on the sign of theta also affect the optimisation result. We analyse this issue via the optimisation of the ground truth square lattice in three cases: θ ≥ 0, θ ≤ 0, and unconstrained sign. When θ ≥ 0, the unique solution is θ 10 ≈ 1 and the generated structure is similar to the target (Fig. 3) . When θ < 0, as shown in Fig. 8 , the only solution is θ 20 ≈ −1, however the structure consists of many 4-, 5-, 6-nearest neighbour particles and is quite irregular. When the sign of θ is not constrained, we are able to obtain a lattice similar to the target in terms of both position and orientation, however the solution is a combination of θ 10 ≈ 0.5 and θ 20 ≈ −0.5. The contribution of θ 10 > 0 to the formation of the square lattice is well confirmed. It facilitates the head to tail alignment of the Y 10 particles, and the side-by-side or parallel alignment. The θ 20 < 0 corresponds to the tri-block particle similar to the one for kagome lattice. The rise of θ 20 < 0 in this case enhances the head-to-tail configuration, and possibly also enhances the parallel alignment because the particle with θ 20 < 0 (i.e. the equatorial parts are attractive) is more energetically favourable. It is interesting that using θ 20 < 0 alone creates a undesired structure, but a combination with another design θ 10 > 0 can create the desired structure. About the choice of sign-constraint, it depends on the situation, for example when the interaction of the red and blue patch is already set, then the constraint is a must. Setting either positive or negative θ often requires less number of iterations because the design parameter space is reduced. The consideration of unconstrained θ may spend more computational cost, but the estimation can lead to richer possible outcomes of patchy particle design.

As mentioned earlier, the absolute value of the coefficient θ refers to the contribution of spherical harmonics to patchy particle design. In the case of square lattice, for example, the estimated result is θ 10 = 1.0 or θ 10 = 1.2 depending on the target. The difference of these values does not matter much, because they both lead to the target structure, and suggests the same particle design.

In the scope of this study, the pool includes five types of patchy particles. It is straightforward to increase the number of parameters. We have increased the number of parameters to 9 by adding four more types of patchy particle Y 30 , Y 31 , Y 32 , Y 33 and performed inverse optimisation for a ground truth square lattice. Regarding the computation cost, the time for each iteration increases by around 30% compared to the 5parameter case. The computational cost scales almost linearly in the number of parameters. This is appealing compare with the conventional grid search, in which the cost scales exponentially. As shown in Fig. 9 , the estimation result of θ is unique and consist of two nonzero parameters θ 10 ≈ 0.65 and θ 30 ≈ 0.25. The estimation is comparable to the 5-parameter case in Fig. 3(a) . The rise of θ 30 is understandable because the Y 30 has similar symmetry to the Y 10 , which promote the head- tail configuration. There is some fluctuation for θ 32 , but such fluctuation does not affect the generated structure. Although the inclusion of many parameters makes the behaviours of the self-assembly and the estimation more complicated, it is expected to reveal more combinations and suggestions on the design of patchy particle. FIG. 9: 9-parameter optimisation for ground truth square lattice. Data are split into 2 groups for better presentation.The generated structure by estimation result, and the estimation of particle patchiness are also given. given particle. Then the orientation of an arbitrary neighbour j , positioned at x j around the centre particle i , denoted aŝ n j (x j |(i , j )) can be interpolated from the density g(n|x) = g n (n, x) g x (x) (D3)

Computational consideration:

(1) About R(n j |n i ): The orientation of particle consists of three local orthonormal basesn (m) , m = 1, 2, 3 in Cartesian coordinates. In principal, when translate and rotate the jth particle so that x i = 0 andn i is on some axis, one can choosen (1) i ,n (2) i ,n (3) i align with Ox, Oy, Oz, respectively; as a result, the position of particle jth after the transformation x j → x ji = x j − x i is in three-dimensional space, which makes the calculation of rotational transformation, Eq. (D1), (D2) andn j (x j |(i , j )) more complex. In this study, since x is in xy plane, we simplify those calculations by choosing the two bases almost lying on the xy plane, then align their xy-projection with the axis Ox, Oy. As a result, the position of particle is always on xy plane, and the rotation n j → R(n j |n i, projected ) is two-dimensional.

(2) For the Gaussian kernel G(x) and G(n), the position of the particle is in polar coordinate x = [0, 2π], while the orientationn is in Cartesian coordinate n = [−1, 1]. The width of Gaussian kernel is chosen so that the full width at half maximum is around 30-50 times smaller than the range of variable.

We choose σ x = 0.1, σ n = 0.02.

(3) Only particles in the 1st shell (i.e. around the first peak of pair distribution function) are considered.

(4) Assign the orientation: As mentioned above, two local bases of particle orientationn (m) ,n (n) are independently estimated, then the 3rd basis is determinedn (l) =n (m) ×n (n) . Since the three bases are required to be orthogonal, singular value decomposition is applied to find the nearest orthogonal matrix for the three bases. During the assignment, the constraintn (m) .n (n) < 0.6 is employed to reduce the deformation of the bases before and after applying SVD. (Note: to find the orthogonal basis, another way is just fixn (m) and rotaten (n) to the nearest orthogonal vector).

(5) The assignment of orientation for target structure can be performed in two ways: global or local. (i) In global assignment, from an initial pair (i, j), the calculation is then propagated to all particles of target. Orientation for jth particle is calculate m( j) times, and we have m( j) configuration for target. Note that we can not take the mean of orientation (i.e. the averagen 2 (x j ) = 1 m( j) ∑ in2 (x j |(i, j)) due to its circular quantity nature. (ii) In local assignment, the calculation is only for neighbours j of particle i.

Assignment of orientations for all particles is not chosen, because the further the neighbours are, the more fluctuating their orientations are.

Self-assembly of block copolymers towards mesoporous materials for energy storage and conversion systems

Self-Assembly of Colloidal Nanocrystals: From Intricate Structures to Functional Materials

Inverse methods for design of soft materials

Directed selfassembly of a colloidal kagome lattice

Designing selfassembling kinetics with differentiable statistical physics models

Colloidal diamond

Computational self-assembly of a one-component icosahedral quasicrystal

Mosaic two-lengthscale quasicrystals

Anisotropy of building blocks and their assembly into complex structures

Predictive Self-Assembly of Polyhedra into Complex Structures

Computational self-assembly of colloidal crystals from Platonic polyhedral sphere clusters

Inverse patchy colloids: Synthesis, modeling and self-organization. Current Opinion in Colloid & Interface Science

Topological defects of dipole patchy particles on a spherical surface

Simple method for the synthesis of inverse patchy colloids

Assembly and Phase Transitions within Colloidal Crystals

Comment: 18 pages, 5 figures, 2 boxes

Probabilistic inverse design for self-assembling materials

Machine learning and data science in soft materials engineering

Optimized Interactions for Targeted Self-Assembly: Application to a Honeycomb Lattice

Designed interaction potentials via inverse methods for selfassembly

Inverse optimization techniques for targeted self-assembly

Design of two-dimensional particle assemblies using isotropic pair interactions with an attractive well

Inverse design of charged colloidal particle interactions for self assembly into specified crystal structures

Engineering entropy for the inverse design of colloidal crystals from hard shapes

Programming patchy particles to form complex periodic structures

Designing Patchy Interactions to Self-Assemble Arbitrary Structures

Learning to grow: Control of material self-assembly using evolutionary reinforcement learning

Role of Surface Charge Density in Nanoparticle-Templated Assembly of Bromovirus Protein Cages

Inverse patchy colloids: From microscopic description to mesoscopic coarse-graining

Computer Simulation of Liquids

Pattern Recognition and Machine Learning. Information Science and Statistics

Aperiodic Crystals

Formation of dodecagonal quasicrystals in twodimensional systems of patchy particles

Activity-Enhanced Self-Assembly of a Colloidal Kagome Lattice

Self-assembly scenarios of patchy colloidal particles in two dimensions

On the stability of Archimedean tilings formed by patchy particles

Minimal Positive Design for Self-Assembly of the Archimedean Tilings

2D Monte Carlo Simulation of Patchy Particles Association and Protein Crystal Polymorph Selection

Two dimensional assembly of triblock Janus particles into crystal phases in the two bond per patch limit

Defects and Geometry in Condensed Matter Physics

The authors acknowledge the support from JSPS KAK-ENHI Grant number JP20K14437 to U.T.L., and JP20K03874 and JP20H05259 to N.Y.

The data that support the findings of this study are available from the corresponding author upon reasonable request.Appendix A: Anisotropic interaction of the patchy particleThe anisotropic interaction is calculated based on the mutual orientation of a pair of particle i and j. Let then (m) i , n (m) j , m = 1, 2, 3 are local bases of particle i and j (Fig.  1) , andr is the unit distance vector between particle center. The anisotropic interaction Ξ lm ∝ {C (l,m)The detail of the isotropic Week-Chandler-Anderson potential u WCA preventing the overlapping of particle, and the Morse potential u M in Eq.(1) is given aswhere r = r i j = r j − r i is the distance vector between particle centres, r = |r|, andr = r/r, ε is the potential well depth, r eq is the Morse potential equilibrium position (r eq = 1.878a), M d , M r is the Morse potential depth and range respectively (M d = 3.761a, M r = 0.213a for short-range potential and M d = 2.294a, M r = a for long-range potential [40] ).Appendix C: Local structure analysisFor each particle j, the number of nearest neighbours n nb is determined by counting the number of particle k satisfying r jk ≤ 2.5a where a is the radius of the particle.The n-fold bond-orientational order parameter [41] ψ n of particle j is calculated by ψ n ( j) = 1 n nb ∑ n nb k=1 e inϕ jk , where ϕ jk is the angle between particle j and its neighbouring particle k. The |ψ n | value characterises the local degree of the regular nfold order around a particle; for example, the hexagonal lattice has |ψ 6 | = 1, kagome lattice also has |ψ 6 | = 1 but the number of nearest neighbours is n nb = 4.

For a given target, information of the target is required. Such information is often the positions of particles, or the density field. In the case of patchy particle, the orientation indicating the patches must be known. In reality the orientation may not be measured. In this case, the orientation must be treated as a hidden variable. Our approach is that the orientation of the target is interpolated from that of the generated structure during each iteration.From the position and orientation of the particles in generated structure, we aim to (i) compute the local orientation field of the generated structure, and (ii) statistically estimate the orientations of particles in targeted structure. From a given generated structure consisting particles of the position {x 1 , x 2 , ...} and orientation {n 1 ,n 2 , ...}, where x ∈ R 2 and n ∈ R 3 in this study, the orientation fieldn(x) is locally estimated. Consider the particle i, the neighbour j can be determined if the distance ||x j − x i || ≤ r nb , and let m(i) be the number of neighbours of particle i. Thenn(x) is taken from the probability distribution P(x j ,n j ) under fixed position and orientation of particle i, i.e. translate and rotate the jth particles so that x i = 0 andn i is on some axis. This is written as x j → x ji = x j − x i andn j → R(n j |n i ). The position of x j is implicitly rotated. Apply kernel density estimator for arbitrary position x in the system, the density of position of j and orientation of j is determined aswhere x ji = x j − x i (after the transformation), N is the number of particle i ∈ [1, N], G(x) is Gaussian kernel (kernel with the shape of a Gaussian curve) defined and normalised as2σ 2 , and . i j is mean over i and j. Here Gaussian kernel is applied for both position x and orientationn in Eq. (D1) and (D2). For the n-dimension data, we apply a one dimensional Gaussian curve sequentially in the n dimensions. The equations (D1) and (D2) includes the density distribution of the neighbours orientations and position around a