# 2.5. Colloids¶

## 2.5.1. Introducing colloids¶

If no relevant key words are present, no colloids will be expected at run time. The simulation will progress in the usual fashion with fluid only.

If colloids are required, the `colloid_init`

key word must be present to allow the code to determine where
colloid information will come from. The options for the
`colloid_init`

key word are summarised below:

```
colloid_init none
# none no colloids [DEFAULT]
# input_one one colloid from input file
# input_two two colloids from input file
# input_three three colloids from input file
# input_random Small number at random
# from_file Read a separate file (including all restarts)
```

For idealised simulations which require 1, 2, or 3 colloids, relevant
initial state information
for each particle can be included in the input file. In principle, most
of the colloid state as defined in the colloid
state structure in `colloid.h`

may be specified. (Some state is
reserved for internal use only and cannot be set from the input file.)
Furthermore, not all the state is relevant in all simulations —
quantities such as charge and wetting parameters may not be required,
in which case they a simply ignored.

A minimal initialisation is shown in the following example:

```
colloid_init input_one
colloid_one_a0 1.25
colloid_one_ah 1.25
colloid_one_r 12.0_12.0_12.0
```

This initialises a single colloid with an input radius \(a_0=1.25\),
and a hydrodynamic radius \(a_h=1.25\); in general both are required,
but they do not have to be equal.
A valid position is required within the extent of the system
\(0.5 < x,y,z < L + 0.5\) as specified by the `size`

key word.
State which is not explicitly defined is initialised to zero.

## 2.5.2. Single colloid initialisation¶

A full list of colloid state-related key words is as follows. All the quantities are floating point numbers unless explicitly stated to be otherwise:

```
# colloid_one_nbonds (integer) number of bonds
# colloid_one_bond1 (integer) index of bond partner 1
# colloid_one_bond2 (integer) index of bond partner 2
# colloid_one_isfixedr colloid has fixed position (integer 0 or 1)
# colloid_one_isfixedv colloid has fixed velocity (integer 0 or 1)
# colloid_one_isfixedw colloid has fixed angular velocity (0 or 1)
# colloid_one_isfixeds colloid has fixed spin (0 or 1)
# colloid_one_shape 'sphere' [the default] or 'ellipsoid'
# colloid_one_bc boundary condition 'bbl' [the default] or 'subgrid'
# colloid_one_active active switch [no]
# colloid_one_magnetic magnetic switch [no]
# colloid_one_a0 input radius
# colloid_one_ah hydrodynamic radius
# colloid_one_al subgrid offset parameter (subgrid only)
# colloid_one_r position vector
# colloid_one_v velocity (vector)
# colloid_one_w angular velocity (vector)
# colloid_one_s spin (unit vector)
# colloid_one_m direction of motion (unit) vector for swimmers
```

Note that for magnetic particles, the appropriate initialisation involves
the spin key word `colloid_one_s`

which relates to the dipole
moment \(\mu \mathbf{s}_i\), while `colloid_one_m`

relates to the
direction of motion vector. Do not confuse the two.
It is possible in principle to have magnetic active particles,
in which case the dipole direction or spin (\(\mathbf{s}_i\)) and the
direction of swimming motion \(\mathbf{m}_i\) are allowed to be distinct.

```
# colloid_one_b1 Squirmer parameter B_1
# colloid_one_b2 Squirmer parameter B_2
# colloid_one_rng (integer) random number generator state
# colloid_one_q0 charge (charge species 0)
# colloid_one_q1 charge (charge species 1)
# colloid_one_epsilon Permeativity
# colloid_one_c Wetting parameter C
# colloid_one_h Wetting parameter H
```

### 2.5.2.1. Example: Single active particle (a squirmer)¶

The following example shows a single active particle with initial swimming direction along the \(x\)-axis.

```
colloid_init input_one
colloid_one_shape sphere
colloid_one_active yes
colloid_one_a0 7.25
colloid_one_ah 7.25
colloid_one_r 32.0_32.0_32.0
colloid_one_v 0.0_0.0_0.0
colloid_one_m 1.0_0.0_0.0
colloid_one_b1 0.05
colloid_one_b2 0.05
```

### 2.5.2.2. Fixed position or velocity¶

It is possible to fix the position or velocity of a colloid via

```
colloid_one_isfixedr 1
colloid_one_isfixedv 1
```

It is also possible to do this on a per co-ordinate direction basis using

```
colloid_one_isfixedrxyz 0_0_1
colloid_one_isfixedvxyz 0_0_1
```

to, for example, fix the \(z\)-poisition and velocity components only.

## 2.5.3. Ellipsoidal particles¶

One may select, e.g.,

```
colloid_one_shape ellipsoid
```

which will allow the use of prolate ellipoids. Additional information is required to define the semi-major axes, and an initial orientation. The semi-major axes \(a, b, c\) are introduced via the key

```
colloid_one_elabc 7.5_2.5_2.5 # a, b, c
```

At the moment that are a number of constraints on this choice. We must have \(a >= b >= c\) for an ordinary ellipsoid, with an additional constraint that \(b = c\) if the particle is active, or requires surface anchoring boundary conditions.

It is not possible to use an ellipse in two dimensions.

### 2.5.3.1. Initial orientation of ellipsoids¶

There are two ways available to specify the initial orientation of an ellipsoid (defined by the lines along the principal axes \(a\) and \(b\)). The first is via a standard \(z-x-z\) set of Euler angles (specified in degrees):

```
colloid_one_euler 90.0_30.0_0.0 # Euler angles (degrees)
```

These angles will be used to compute the initial orientation internally.

The second way to be define two vectors with align along the first and second principle axis (\(a\) and \(b\)). This is done via

```
colloid_one_elev1 1.0_1.0_1.0 # First vector
colloid_one_elev2 1.0_0.0_0.0 # Second vector
```

The vectors do not have to be unit vectors; they do not even have to be at right angles (the second will be orthogonalisated against the first). However, they must no be parallel (linearly dependent). Again, the related orientation (including the third axis) will be computed internally.

### 2.5.3.2. Initialisation of orientation in colloid state files¶

The internal representation of current ellipsoid orientation is using a unit quaternion. This is what appears in the colloid state file used for input and output. If one wants to generate an appropriate input file, the relevant unit quaternion must be computed.

## 2.5.4. Many colloid initialisation at random¶

For suspensions with more than few colloids, but still at relatively low volume fraction (10–20% by volume), it is possible to request initialisation at random positions.

The additional key word value pair `colloid_random_no`

determines the total number of particles to be placed in
the system. To prevent particles being initialised very
close together, which can cause problems in the first few
time steps if strong potential interactions are present,
a ```
grace'' distance or minimum surface-surface separation
may also be specified (``colloid_random_dh
```

).

The following example asks for 100 colloids to be initialised at random positions, with a minimum separation of 0.5 lattice spacing.

```
colloid_init input_random
colloid_random_no 100 # Total number of colloids
colloid_random_dh 0.5 # ``Grace'' distance
colloid_random_a0 2.30
colloid_random_ah 2.40
```

An input radius and hydrodynamic radius must be provided: these
are the same for all colloids.
If specific initialisations of the colloid state (excepting the
position) other than the radii are wanted, values should be provided
as for the single particle case in the preceding section, but using
`colloid_random_a0`

in place of `colloid_one_a0`

and so on.

The code will try to initialise the requested number in the current system size, but only makes a finite number of attempts to place particles at random with no overlaps. (The initialisation will also take into account the presence of any solid walls, using the same grace distance.) If the the number of particles is too large, the code will halt with a message to that effect.

In general, colloid information for a arbitrary configuration with many particles should be read from a pre-prepared file. See the section on File I/O for further information on reading files.

## 2.5.5. Colloid interactions¶

Note that two-body pair-potential interactions are defined uniformly for all colloids in a simulation. The same is true for lubrication corrections. There are a number of constraints related to the computation of interactions discussed below.

### 2.5.5.1. Boundary-colloid lubrication correction¶

Lubrication corrections (here the normal force) between a flat wall are required to prevent overlap between colloid and the wall. The cutoff distance is set via the key word value pair

```
boundary_lubrication_rcnormal 0.5 # normal cut-off distance
# Note this is 'rcnormal'
```

It is recommended that this is used in all cases involving walls. A reasonable value for the cut off is in the range \(0.1 < h_c < 0.5\) in lattice units, and should be calibrated for particle hydrodynamic radius and fluid viscosity \(\eta\) if exact results are important.

The form of the lubrication correction for a colloid of hydrodynamic radius \(a_h\) and velocity \(U_\alpha\) is based on the analytical expression for the lubrication force between a place wall and a sphere:

where \({\hat n}_\alpha\) is the unit normal between the wall and the centre of the colloid. The correction is zero for \(h > h_c\). The surface-surface separation is \(h\) and the cut off is \(h_c\) as illustrated in the following diagram on the left. In this illustration, the wall position (full line) is at the default value of \(x = 0.5`for the lower end of the system in the :math:`x\)-direction.

In some situations it may be useful to prevent colloids approaching the wall to within some fixed separation. It is possible to provide a uniform offset the apparent position of the wall by a fixed amount \(\delta h\). This is illustrated on the right above. The lubrication correction will then diverge at the new wall position (dotted line). This can be useful to maintain clear fluid lattice sites between wall and colloid.

The relevant input key is

```
boundary_lubrication_dhnormal 0.2 # normal offset distance (+ve)
# [Default: 0]
```

For both the cut-off, and the offset distance, there is a single input value which takes effect in all three co-ordinate directions if a wall is present in the corresponding direction. In practice, the lubrication correction between wall and colloid should be very robust.

### 2.5.5.2. Boundary-colloid soft sphere potential¶

In some circumstances it may be desirable to use a conservative potential at a boundary wall in place of the lubrication correction. In this case a cut-and-shifted soft sphere potential is available. Foe example:

```
wall_ss_cut_on yes # Switch
wall_ss_cut_epsilon 0.001 # Energy scale
wall_ss_cut_sigma 0.1 # Length scale
wall_ss_cut_nu 2.0 # Exponent
wall_ss_cut_hc 0.5 # wall-surface cut off
```

Both the exponent and the wall-surface cut off should be positive. The potential will take effect at boundary walls in all directions.

### 2.5.5.3. Colloid-colloid lubrication corrections¶

The key words to activate the calculation of lubrication corrections are:

```
lubrication_on 1
lubrication_normal_cutoff 0.5
lubrication_tangential_cutoff 0.05
```

### 2.5.5.4. Soft-sphere potential¶

A cut-and-shifted soft-sphere potential of the form \(v \sim \epsilon (\sigma/r)^\nu\) is available. Some trial-and-error with the parameters may be required in any given situation to ensure simulation stability in the long run. The following gives an example of the relevant input key words:

```
soft_sphere_on 1 # integer 0/1 for off/on
soft_sphere_epsilon 0.0004 # energy units
soft_sphere_sigma 1.0 # a length
soft_sphere_nu 1.0 # exponent is positive
soft_sphere_cutoff 2.25 # a surface-surface separation
```

### 2.5.5.5. Soft-sphere potential (type-specific)¶

This potential is of the same form as the basic cut-and-shifted
soft-sphere potential
described above, but allows different parameters to be specified for
colloids with different *interaction type*. The interaction type is
an integer specified by the appropriate element of the colloid
structure, e.g., via input

```
colloid_one_interact_type 0
...
colloid_two_interact_type 1
```

specifying two different types (0 and 1). The first type must have index 0. Interactions between different pairs then all have the form \(v_{ij} \sim \epsilon_{ij} (\sigma_{ij}/r)^{\nu_{ij}}\).

The type specific pair interaction is then introduced via

```
pair_ss_cut_ij yes
pair_ss_cut_ij_ntypes 2
```

the second key value pair giving the number of types expected. The parameters then form a symmetric matrix, for which we specific the upper triangle as a flattened vector. In the case of two types, there are three independent parameters, e.g.,

```
pair_ss_cut_ij_epsilon 0.2_0.1_0.05 # epsilon_00, _01, _11 in order
```

where we specify \(\epsilon_{00}, \epsilon_{01}\) and \(\epsilon_{11}\), being the interaction energies for interactions bewtween pairs of type (0,0), (0,1), and (1,1) respectively. The value \(\epsilon_{10}\) is set to be the same as \(\epsilon_{01}\) internally. A full set of key value pairs might be

```
pair_ss_cut_ij yes # Switch on
pair_ss_cut_ij_ntypes 2 # Number of types n
pair_ss_cut_ij_epsilon 0.0_0.1_0.0 # n(n+1)/2 epsilon parameters
pair_ss_cut_ij_sigma 0.0_2.0_0.0 # n(n+1)/2 sigma parameters
pair_ss_cut_ij_nu 1.0_1.0_3.0 # n(n+1)/2 nu exponents
pair_ss_cut_ij_hc 0.1_0.4_0.1 # n(n+1)/2 surface-surface cut offs
```

The user must ensure all colloids have appropriate interaction types, i.e., the interaction type does not exceed 1 in this case.

### 2.5.5.6. Lennard-Jones potential¶

The Lennard-Jones potential is controlled by the following key words:

```
lennard_jones_on 1 # integer 0/1 off/on
lj_epsilon 0.1 # energy units
lj_sigma 2.6 # potential length scale
lj_cutoff 8.0 # a centre-centre separation
```

### 2.5.5.7. Yukawa potential¶

A cut-and-shifted Yukawa potential of the form \(v \sim \epsilon \exp(-\kappa r)/r\) is available using the following key word value pairs:

```
yukawa_on 1 # integer 0/1 off/on
yukawa_epsilon 1.330 # energy units
yukawa_kappa 0.725 # an inverse length
yukawa_cutoff 16.0 # a centre-centre cutoff
```

### 2.5.5.8. Dipole-dipole interaction and the Ewald sum¶

The Ewald sum is completely specified in the input file by the uniform dipole strength $mu$ and the real-space cut off \(r_c\).

```
ewald_sum 1 # integer 0/1 off/on
ewald_mu 0.285 # dipole strength mu
ewald_rc 16.0 # real space cut off
```

If short range interactions are required, particle information is stored in a cell list, which allows efficient computation of the potentially \(N^2\) interactions present. This gives rise to a constraint that the width of the cells must be large enough that all relevant interactions are included. This generally means that the cells must be at least \(2a_h + h_c\) where \(h_c\) is the largest relevant cut off distance.

The requirement for at least two cells per local domain in parallel means that there is a associated minimum local domain size. This is computed at run time on the basis of the input. If the local domain is too small, the code will terminate with an error message. The local domain size should be increased.

### 2.5.5.9. External forces¶

The following example requests a uniform body force in the negative \(z\)-direction on all particles.

```
colloid_gravity 0.0_0.0_-0.001 # vector
```

The counterbalancing body force on the fluid which enforces the constraint of momentum conservation for the system as a whole is computed automatically by the code at each time step.

Note: in a real system, a gravitational force on a colloid is related to buoyancy \(F \propto \Delta\rho g\), where the density difference is that between the colloid and the surrounding fluid, and \(g\) is an acceleration. In a system where there is no density contrast, as we have here (typically), the “gravity” is the product \(\Delta\rho g\). Formally, this may be viewed as the limit that \(\Delta\rho \rightarrow 0\), combined with the limit \(g \rightarrow \infty\), but the limit of the product is finite.

As a separate alternative, one can specify

```
colloid_buoyancy 0.0_0.0_-0.001 # vector
```

which introduces a force proportional to the volume of each individual
colloid. (The `colloid_gravity`

option is a constant force, independent
of colloid size.) One cannot
specify both a buoyancy and a gravity force at the same time. The
appropriate counterbalancing force on the fluid is again computed
internally.

One should not specify both a gravitational force (or buoyancy) at the same time a a body force on the fluid (see below).

Note that if the system has a “bottom”, e.g., a plane wall normal to the direction of the gravitational force, the counterbalancing force on the fluid is not really required. However, it is always present in the current implementation irrespective of any walls.

### 2.5.5.10. Body force on the fluid when colloids are present¶

It is possible to impose an external body force on the fluid as described in Constant body forces to provide, effectively, a pressure gradient. If colloids are present, a contribution to the force on the colloid proportional to the colloid’s discrete volume is applied. This gives the same total momumtum input as if the colloid were replaced by fluid.

It is not possible to have both a both an external gravitational force on
the colloid (as via `colloid_gravity`

) and a body force on the fluid at
the same time.

## 2.5.6. Liquid crystal anchoring at colloid surfaces¶

The preferred orientation of the liquid crystal director at the surface of a colloid can be of one of two different types:

```
lc_anchoring_coll normal
lc_anchoring_coll planar
```

For both cases, the chosen condition appliess to all colloids in the system in the same way.

The liquid crystal anchoring boundary condition is implemented via the calculation of the tensor order parameter gradients at the surface of the colloid. We assume there is a surface free energy density (per unit area)

where \(Q_{\alpha\beta}\) us the adjacent fluid order parameter, and \(Q^0_{\alpha\beta}\) is some preferred otder parameter determined by the type of anchoring.

The boundary condition is derived from the Euler-Lagrange equation, and contains gradient terms in the bulk free energy density and the surface free energy density \(f\) and \(f_s\), along with the outward unit normal at the surface \(\hat{n}_\gamma\):

A suitable gradient computation must be selected (see below).

### 2.5.6.1. Normal (or homoetropic) anchoring at colloid surfaces¶

If the preferred orientation of order at the colloid surface is normal to the surface (also referred to as homoetropic anchoring), the surface free energy density per unit area may be written

Here, \(w_1\) is a parameter, and \(Q_{\alpha\beta}\) is the local fluid order parameter. The preferred orientation is based on the outward normal \(\hat{n}_\gamma\) at the surface (in the direction from the centre of the colloid to the relevant fluid site at the surface). The order parameter tensor

The value of \(A\) corresponds to that which minimises the bulk free energy (depending on the temperature \(\gamma\)). The full boundary condition for the order parameter gradient at the colloid surface then contains the term

The relevant input key/value pairs for normal colloid anchoring are:

```
lc_coll_anchoring normal # anchoring type
lc_coll_anchoring_w1 0.002 # free energy parameter w_1
```

It is often appropriate to set the value of the surface free energy parameter in the context of the bulk elastic constant, e.g., by considering the dimensionless group \(w_1/\kappa a\), where \(a\) is the radius of the colloid.

As an example, an ordinary nematic is initialised with the director along the \(x\)-direction, and a single stationary sphereical colloid with normal anchoring introduced to the system. The nematic is allowed to relax with no hydrodynamics. A section of the resulting director field is shown for a “weak” anchoring case (left), and a “strong” anchoring case (right). The \(x\)-direction is in the horizontal. In the strong anchoring case, it can be seen that a defect in the orientational order has appeared above and below the surface of the colloid. In three dimensions, this defect is present all around the circumference of the colloid and forms a “Saturn ring”.

While the nominal radius of the colloid is indicated by the circle, it should be remembered that the discrete shape is block-like, as indicated by the presence of director in individual grid cells.

### 2.5.6.2. Planar (or degenerate) anchoring at colloid surfaces¶

For planar anchoring, the preferred orientation is in the local tangent plane at the surface: this is a degenerate case as any orientation in the plane is energetically equivalent. An appropriate boundary condition is described by Fournier and Galatola [FournierGalatola2005], which we write as

To compute this term we take the local fluid order parameter \(Q_{\alpha\beta}\), form the quantity

which is then projected onto the tangent plane via \(\tilde{Q}^\perp_{\alpha\beta} = P_{\alpha\gamma} \tilde{Q}_{\gamma\sigma} P_{\sigma\beta}\) with the local surface normal entering through \(P_{\alpha\beta} = \delta_{\alpha\beta} - \hat{n}_\alpha \hat{n}_\beta\). (Again, the normal at the colloid surface is based on the displacement from the colloid centre.) The full boundary condition arising from the surface free energy contains the terms

The term \(S_0 = 3A/2\), with amplitude \(A\) as described above for normal anchoring.

Relevant input parameters are:

```
lc_coll_anchoring planar
lc_coll_anchoring_w1 0.01 # both w1 and w2 must be present
lc_coll_anchoring_w2 0.005
```

Following the example for normal anchoring, the illustration below shows the result for planar anchoring with \(w_1 = w_2\) (weak and strong cases are left and right, respectively). In the strong case one can identify a pair of point defects at either side of the colloid, usually referred to as “boojums”. Note that the director has actually rotated into the third dimension at these points and so appears forshortened.

### 2.5.6.3. Anchoring when more than one colloid is present¶

If more than one colloid is present in the system, then the surface normal can become poorly defined if adjacent lattice sites are occupied by different colloids. If such cases the anchoring term is set to zero in the boundary condition. The advice here is to prevent close approaches between colloids by means of, e.g., a soft-sphere potential. This should ensure that fluid sites are always present in the gap.