hydro_util

Short description

@author: Moritz F P Becker

pyrid.molecules.hydro_util.A(r)

Returns matrix A which turns the cross product rxw into a dot product A.w.

Parameters

rfloat

Bead position

Notes

We can turn the cross product \(r \times \omega\) with \(r = [x,y,z]\) into a dot product \(A \cdot \omega\), where

\[\begin{split}A = \begin{pmatrix} 0 & -z & y\\ z & 0 & -x\\ -y & x & 0 \end{pmatrix}.\end{split}\]

pyrid.molecules.hydro_util.P(rij)

Calculates the normalized tensor/outer product \(\frac{\boldsymbol{r}_{ij} \otimes \boldsymbol{r}_{ij}}{r_{ij}}\) for the hydrodynamic interaction tensor (Oseen tensor).

Parameters

rijfloat64[3]

Distance between partcile i and j

Returns

float[3,3]

Normalized outer product of r_ij

pyrid.molecules.hydro_util.calc_CoD(D_rr, D_tr)

Calculates the center of diffusion of a rigid bead molecule.

Parameters

D_rrfloat64[3,3]

Rotational diffusion tensor.

D_trfloat64[3,3]

Translation-rotation coupling

Returns

float64[3]

Center of diffusion.

Notes

The diffusion tensor \([[D^{tt}, D^{tr}],[D^{rt},D^{rr}]]\) depends on the choice of the bead model origin (body frame). For symmetry reason, one should use the so-called center of Diffusion, which can be calculated from a diffusion tensor referring to an arbitrary origin [15], [2] :

(7)\[\begin{split}\begin{align*} \boldsymbol{r}_{OD} = & \begin{pmatrix} x_{OD} \\ y_{OD}\\ z_{OD} \end{pmatrix} \\ = & \begin{pmatrix} D_{rr}^{yy}+D_{rr}^{zz} & -D_{rr}^{xy} & -D_{rr}^{xz}\\ -D_{rr}^{xy} & D_{rr}^{xx}+D_{rr}^{zz} & -D_{rr}^{yz}\\ -D_{rr}^{xz} & -D_{rr}^{yz} & D_{rr}^{yy}+D_{rr}^{xx} \end{pmatrix}^{-1} \begin{pmatrix} D_{tr}^{zy}-D_{tr}^{yz}\\ D_{tr}^{xz}-D_{tr}^{zx}\\ D_{tr}^{yx}-D_{tr}^{xy} \end{pmatrix} \end{align*}\end{split}\]

After calculating the the center of diffusion, we simply set the rigid body origin to \(r_{OD}\) and recalculate the diffusion tensor.

pyrid.molecules.hydro_util.calc_D(Xi_tt, Xi_rt, Xi_tr, Xi_rr, kB, Temp)

Caclulates the diffusion tensor from the supermatrix inverse of the friction tensor.

Parameters

parameter_1dtype

Some Information

parameter_2dtype

Some Information

Returns

tuple(float64[3,3], float64[3,3], float64[3,3], float64[3,3])

Translation, rotational and rotation-translation coupling diffusion tensors D_tt, D_tr, D_rt, D_rr.

Notes

The diffusion tensor can be calculated from the rigid body’s friction tensor via [2] :

(8)\[\begin{split}\begin{pmatrix} \boldsymbol{D}^{tt} & \boldsymbol{D}^{tr,T} \\ \boldsymbol{D}^{rt} & \boldsymbol{D}^{rr} \\ \end{pmatrix} = k_B T \begin{pmatrix} \boldsymbol{\Xi}^{tt} & \boldsymbol{\Xi}^{tr,T} \\ \boldsymbol{\Xi}^{rt} & \boldsymbol{\Xi}^{rr} \\ \end{pmatrix}^{-1} \end{split}\]

pyrid.molecules.hydro_util.calc_Xi(zeta_tt, zeta_rt, zeta_tr, zeta_rr, pos, radii)

Calculates the friction tensors (3,3) of a ridig bead molecule from the friction super matricies (3N,3N) of its individiual beads.

Parameters

zeta_ttfloat64[3N,3N]

Translational friction matrix.

zeta_rtfloat64[3N,3N]

Rotation-translation coupling matrix.

zeta_trfloat64[3N,3N]

Translation-rotation coupling matrix.

zeta_rrfloat64[3N,3N]

Rotational friction matrix.

my_moleculeobj

Molceule class instance containing position data of each bead.

Returns

tuple(float[3,3], float[3,3], float[3,3], float[3,3])

Returns the friction tensors Xi_tt, Xi_rt, Xi_tr, Xi_rr of the molecule.

Raises

NotImplementedError (just an example)

Brief explanation of why/when this exception is raised

Notes

Here, we closely follow [16]. To get an expression for the friction tensor of a rigid bead molecule, we start by considering a sytsem of \(N\) free spherical particles in a fluid with viscosity \(\eta_0\). Each sphere lateraly moves at some velocity :math`u_i` and rotates with some angular velocity \(\omega_i\). The spheres will experience a frictional force and torque \(F_i, T_i\). In the noninertial regime (Stokes), the relationship between the force/torque and the velocities are linear:

(9)\[\boldsymbol{F}_i = \sum_{j=1}^N \zeta_{ij}^{tt} \cdot \boldsymbol{u}_j + \zeta_{ij}^{tr} \cdot \boldsymbol{\omega}_j\]

(10)\[\boldsymbol{T}_i = \sum_{j=1}^N \zeta_{ij}^{rt} \cdot \boldsymbol{u}_j + \zeta_{ij}^{rr} \cdot \boldsymbol{\omega}_j .\]

The \(\zeta_{ij}\) are the (3x3) friction matrices, connecting the amount of friction a particle i expiriences due to the presence of particle j moving through the fluid with velocities \(u_j, \omega_j\). We may rewrite this in a matrix representation as:

(11)\[\begin{split}\begin{pmatrix} F \\ T \\ \end{pmatrix} = \begin{pmatrix} \zeta^{tt} & \zeta^{tr} \\ \zeta^{rt} & \zeta^{rr} \\ \end{pmatrix} \begin{pmatrix} U \\ W \\ \end{pmatrix},\end{split}\]

where \(F = (\boldsymbol{F}_1, ..., \boldsymbol{F}_N)^T\), \(T = (\boldsymbol{T}_1, ..., \boldsymbol{T}_N)^T\) and \(U = (\boldsymbol{u}_1, ..., \boldsymbol{u}_N)^T\), \(W = (\boldsymbol{\omega}_1, ..., \boldsymbol{\omega}_N)^T\). Here \(\zeta\) are of dimension (3Nx3N), forming the friction supermatrix of dimension (6N,6N). The inverted friction supermatrix is the mobility supermatrix (for inmversion of supermatrices also see supermatrix_inverse()).

Next, we consider not a system of N free beads but a rigid bead model, i.e. the beads are rigidly connected. Thereby, all beads move together with some translational velocity \(u_{O}\). Let the body’s frame of reference lie at the center of diffusion of the bead model \(\boldsymbol{r}_O\) and let \(\omega\) be the angular velocity of the rigid bead model. Then, in addition to the translational velocity of the molecule’s center, each bead experiences a translation velocity due to the rotation \(\boldsymbol{\omega} \times \boldsymbol{r}_i\), where \(\boldsymbol{r}_i\) is the position vector from the moclules origin \(\boldsymbol{r}_O\) (in the body frame of reference). Thereby, the total velocity is:

(12)\[\boldsymbol{u}_i = \boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_i \]

Thereby, the force that a single bead experiences due to the movement of all the other beads is:

(13)\[\boldsymbol{F}_i = \sum_{j=1}^N \zeta_{ij}^{tt} \cdot (\boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_j) + \zeta_{ij}^{tr} \cdot \boldsymbol{\omega},\]

and the torque that single bead experiences due to the movement of all the other beads is:

(14)\[\boldsymbol{T}_{P,i} = \sum_{j=1}^N \zeta_{ij}^{rt} \cdot (\boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_j) + \zeta_{ij}^{rr} \cdot \boldsymbol{\omega} . \]

From these expressions we simply get the total force acting at the rigid body origin by summation over all beads:

(15)\[\boldsymbol{F} = \sum_{i=1}^N \sum_{j=1}^N \zeta_{ij}^{tt} \cdot (\boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_j) + \zeta_{ij}^{tr} \cdot \boldsymbol{\omega}\]

For the total torque, however, we get an extra term. \(\boldsymbol{T}_{P,i}\) is only the torque acting on bead i relative to it’s center, i.e. the center of the sphere. Thereby, this only describes the amount of rotation bead i would experience around its center due to the movement of all the other beads. However, the force \(\boldsymbol{F}_{i}\) acting on bead i due to the movement of the other beads also results in a torque with which bead i acts on the rigid bead models center \(r_O\):

(16)\[\boldsymbol{r}_i \times \boldsymbol{F}_i = \boldsymbol{r}_i \times \Big( \sum_j^N \zeta_{ij}^{tt} (\boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_j) + \zeta_{ij}^{tr} \omega \Big)\]

Thereby, the total torque acting on the rigid bead model’s origin is:

(17)\[\boldsymbol{T}_O = \sum_i^N \boldsymbol{T}_{P,i} + \boldsymbol{r}_i \times \boldsymbol{F}_i = \sum_{i=1}^N \sum_{j=1}^N \zeta_{ij}^{rt} \cdot (\boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_j) + \zeta_{ij}^{rr} \cdot \boldsymbol{\omega} + \boldsymbol{r}_i \times \Big( \zeta_{ij}^{tt} (\boldsymbol{u}_O + \boldsymbol{\omega} \times \boldsymbol{r}_j) + \zeta_{ij}^{tr} \omega \Big). \]

In principle, we are done now, however, we would like to transform this into a more ‘general’ expression that we can write in a simple matrix form. For this, we use a little trick to get rid of the cross product, by turning \(\omega \times r\) into the dot product \(- A \cdot \omega\) (note: the sign changed, because of the anticommutativity of the cross product). After some rearranging, we end up with:

(18)\[\boldsymbol{F} = \Big( \sum_{i=1}^N \sum_{j=1}^N \zeta_{ij}^{tt} \Big) \cdot \boldsymbol{u}_O + \Big( \sum_{i=1}^N \sum_{j=1}^N - \zeta_{ij}^{tt} \cdot \boldsymbol{A}_j + \zeta_{ij}^{tr} \Big) \cdot \boldsymbol{\omega}\]

(19)\[\boldsymbol{T} = \Big( \sum_{i=1}^N \sum_{j=1}^N \zeta_{ij}^{rt} + A_i \zeta_{ij}^{tt} \Big) \cdot \boldsymbol{u}_O + \Big( \sum_{i=1}^N \sum_{j=1}^N \zeta_{ij}^{rt} \cdot \boldsymbol{A}_j + \zeta_{ij}^{rr} - A_i \zeta_{ij}^{tt} A_j + A_i \zeta_{ij}^{tr} \Big) \cdot \boldsymbol{\omega}. \]

If we now want write this in matrix form, similar to the free bead example from above:

(20)\[\begin{split}\begin{pmatrix} \boldsymbol{F} \\ \boldsymbol{T}_O \\ \end{pmatrix} = \begin{pmatrix} \Xi^{tt} & \Xi^{tr} \\ \Xi^{rt} & \Xi^{rr} \\ \end{pmatrix} \begin{pmatrix} \boldsymbol{u}_O \\ \boldsymbol{\omega} \\ \end{pmatrix}, \end{split}\]

Where we call \(\Xi\) the friction tensor of the rigid bead molecule [16] :

(21)\[\begin{split}\begin{align*} &\Xi^{tt} = \sum_{i=1}^N \sum_{j=1}^N \zeta_{ij}^{tt} \\ &\Xi_{O}^{tr} = \sum_{i=1}^N \sum_{j=1}^N ( -\zeta_{ij}^{tt} \cdot \boldsymbol{A}_j + \zeta_{ij}^{tr} ) \\ &\Xi_{O}^{rt} = \sum_{i=1}^N \sum_{j=1}^N ( \boldsymbol{A}_j \cdot \zeta_{ij}^{tt} + \zeta_{ij}^{rt} ) \\ &\Xi_{O}^{rr} = \sum_{i=1}^N \sum_{j=1}^N ( \zeta_{ij}^{rr} - \zeta_{ij}^{rt} \cdot \boldsymbol{A}_j + \boldsymbol{A}_i \cdot \zeta_{ij}^{tr} - \boldsymbol{A}_i \cdot \zeta_{ij}^{tt} \boldsymbol{A}_j) \end{align*}\end{split}\]

The only thing left to do now, is calculating \(\zeta\), which we get from the inverse of the mobility supermatrix. The mobility supermatrix is calculated in calc_mu(my_molecule, eta_0).

pyrid.molecules.hydro_util.calc_mu(pos, radii, eta_0)

Calculates the transtalional, rotational and translation-rotation coupling parts of the mobility super matrix \(\mu_{tt}, \mu_{rr}, \mu_{tr} = \mu_{rt}\).

Parameters

my_moleculeobj

Instance of molecule containing the position data of each bead.

eta_0float

Viscosity

Returns

tuple(float[N,N,3,3], float[N,N,3,3], float[N,N,3,3], float[N,N,3,3])

mu_tt, mu_rt, mu_tr, mu_rr

Notes

The mobility matrices are directly related to the modified Oseen tensor. The Oseen tensor has first been introduced by Oseen in 1927 (For reference also see [17]). The Oseen tensor comes up when solving the Stokes equations, which are a linearization of the Navier-Stokes equations. More precisely, it comes up when solving for the flow velocity field in case of a force acting on a point particle (\(f(r) = f_0 \delta(r-r_p)\)) which is immersed in a viscous liquid. In this case, we can write the solution to the Stokes equation as (due to its linearity, any solution to the Stokes equation has to be a linear transformation):

\[v(r) = T(r-r_p) \cdot F\]

T is called the hyrodynamic interaction tensor, Oseen tensor or Green’s function of the Stoke’s equations. The above solution is also called Stokeslet:

\[\boldsymbol{T}(r) = \frac{1}{8 \pi \eta r} \cdot \Big(\boldsymbol{I}+\frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2} \Big),\]

where \(\eta\) is the viscosity and \(\otimes\) is the outer product. In prosa, \(T\) gives the fluid flow velocty at some point \(r\), given a force acting at another point \(r_p\). Kirkwood and Riseman calculated the translational mobility tensor of a rigid bead molecule using the Oseen tensor to describe the hydrodynamic interaction between the beads, and by assigning each bead its friction coefficient \(\zeta_i = 6 \pi \eta_0 a_i\):

\[\begin{split}\begin{align*} \mu_{ij}^{tt} = & \delta_{ij}(6 \pi \eta_0 a_i)^{-1} \boldsymbol{I} \\ & + (1-\delta_ij)(8 \pi \eta_0 r_{ij})^{-1} \\ & \Big(\boldsymbol{I}+\frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2} \Big) \end{align*}\end{split}\]

This solution is fairly intuitive. The first term is just the mobility of a single particle with radius \(a_i\). The second term is just the Oseen tensor. Recalling that the mobility \(\mu\) is defined as the ratio of a particles drift velocity and the applied applied force, the interpretation of the Oseen tensor as representing the interaction part of the bead mobility matrix feels natural (recalling its origin (see above)). However, we also instantly see that something is missing since the Oseen tensor only considers the distance between the bead centers but neglects their volume/radius \(a_i\). Fortunately, Torre and Bloomfield [18] established a correction to the Oseen tensor for nonidentical spheres (also see [19]):

(22)\[\boldsymbol{T}_{ij} = \frac{1}{8 \pi \eta r} \cdot \Big(\boldsymbol{I}+\frac{\boldsymbol{r}_{ij} \otimes \boldsymbol{r}_{ij}}{r_{ij}^2} + \frac{\sigma_i + \sigma_j}{r_{ij}^2} \Big( \frac{1}{3} \boldsymbol{I} - \frac{\boldsymbol{r}_{ij} \otimes \boldsymbol{r}_{ij}}{r_{ij}^2} \Big) \Big),\]

By that, the friction tensor reads:

(23)\[\begin{split}\begin{align*} \mu^{tt}_{ij} = & \delta_{ij} (6 \pi \eta_0 a_i)^{-1} \boldsymbol{I} + (1-\delta_{ij})(8 \pi \eta_0 r_{ij}^{-1})(\boldsymbol{I}+\boldsymbol{P}_{ij}) \\ & + (8 \pi \eta_0 r_{ij}^{-3})(a_i^2+a_j^2)(\boldsymbol{I}-3 \boldsymbol{P}_{ij}), \end{align*}\end{split}\]

where \(\boldsymbol{P}_{ij} = \Big(\boldsymbol{I}+\frac{\boldsymbol{r} \otimes \boldsymbol{r}}{r^2} \Big)\). The mobility tensor for rotation, not correcting for the beads volume, reads [16].

(24)\[\begin{split}\begin{align*} \mu^{rr}_{ij} = & \delta_{ij} (8 \pi \eta_0 a_i^3)^{-1} \boldsymbol{I} \\ & + (1 - \delta_{ij})(16 \pi \eta_0 r^3_{ij})^{-1} (3 \boldsymbol{P}_{ij} - \boldsymbol{I}) \end{align*} \end{split}\]

Here, again, the first term is just the rotational mobility of the single bead and the second term accounts for the hydrodynamic interaction. In this formulation, there is still a correction for the volume missing. This correction consists of adding \(6 \eta_0 V_m \boldsymbol{I}\) to the diagonal components of the rotational friction tensor \(\Xi^{rr}_O\), where \(V_m\) is the volume of the bead model (sum over all bead volumes) [20], [16].

And, at last, for rotation-translation coupling, we have [16]:

(25)\[\mu^{rt}_{ij} = (1-\delta_{ij}) (8 \pi \eta_0 r_{ij}^2)^{-1} \boldsymbol{\epsilon}\boldsymbol{\hat{r}}_{ij} \]

pyrid.molecules.hydro_util.calc_zeta(mu_tt, mu_tr, mu_rt, mu_rr)

Calculates the friction tensors \(\zeta^{tt}, \zeta^{rr}, \zeta^{tr}, \zeta^{rt}\) from the inverse of the mobility supermatrix \([[\mu^{tt}, \mu^{rr}], [\mu^{tr}, \mu^{rt}]]\).

Parameters

mu_ttfloat64[3N,3N]

Translational mobility tensor.

mu_trfloat64[3N,3N]

Translation-rotation coupling.

mu_rtfloat64[3N,3N]

Rotation-translation coupling.

mu_rrfloat64[3N,3N]

Rotational mobility tensor.

Returns

tuple(float64[3N,3N], float64[3N,3N], float64[3N,3N], float64[3N,3N])

The friction tensors \(\zeta^{tt}, \zeta^{rr}, \zeta^{tr}, \zeta^{rt}\).

Raises

NotImplementedError (just an example)

Brief explanation of why/when this exception is raised

pyrid.molecules.hydro_util.center_of_mass(pos, radii)

Returns the center of mass of the molecule assuming equal mass density for all paricles with radius a_i.

Parameters

my_moleculeobj

Instance of molecule containing the position data of each bead.

eta_0float

Viscosity

Returns

float64[3]

Center of mass

Notes

The center of mass is defined by:

(26)\[\boldsymbol{R} = \frac{1}{M} \sum_{i=1}^n m_i \boldsymbol{r_i},\]

where \(M\) is the total mass of the molecule.

pyrid.molecules.hydro_util.diffusion_tensor(Simulation, molecule_name, return_CoD=False, return_CoM=False, return_coupling=False)

Calculates the diffusion tensor, accounting for the center of diffusion of the rigid bead molecule. The origin of the molecule is automatically updated to the center of diffusion \(r_{OD}\)!

Parameters

my_moleculeobj

Instance of molecule containing the position data of each bead.

Simulationobj

Instance of the Simulation class

Returns

tuple(float64[3,3], float64[3,3], float64[3,3], float64[3,3], float64)

Translation, rotational and rotation-translation coupling diffusion tensors D_tt, D_rr, D_tr, D_rt and the center of diffusion r_CoD, center of mass r_CoM.

pyrid.molecules.hydro_util.diffusion_tensor_off_center(pos, radii, eta_0, Temp, kB)

Calculates the diffusion tensor (off the center of diffusion) from the moelcule structure, the fluid viscosity and temperature.

Parameters

my_moleculeobj

Instance of molecule containing the position data of each bead.

eta_0float

Viscosity

Tempfloat

Temperature

kBfloat

Boltzmann constant

Returns

tuple(float64[3,3], float64[3,3], float64[3,3], float64[3,3])

Translation, rotational and rotation-translation coupling diffusion tensors D_tt, D_tr, D_rt, D_rr.

pyrid.molecules.hydro_util.levi_civita(rij)

Returns the product of the Levi-Civita tensor with vector rij.

Parameters

rijfloat64[3]

Distance between partcile i and j

Returns

float[3,3]

Levi-Civita tensor

Notes

Given the distance vector \(r_{ij} = [x_{ij},y_{ij},z_{ij}]\) the product \(\epsilon \cdot r_{ij}\) is returned, where \(\epsilon\) is the Levi-Civita tensor:

\[\begin{split}\epsilon \cdot r_{ij} = \begin{pmatrix} 0 & z_{ij} & -y_{ij}\\ -z_{ij} & 0 & x_{ij}\\ y_{ij} & -x_{ij} & 0 \end{pmatrix}.\end{split}\]

pyrid.molecules.hydro_util.supermatrix_inverse(M1, M2, M3, M4)

Calculates the inverse of a 2x2 supermatrix.

Parameters

M1float64[3N,3N]

submatrix at (0,0)

M2float64[3N,3N]

submatrix at (0,1)

M3float64[3N,3N]

submatrix at (1,0)

M4float64[3N,3N]

submatrix at (1,1)

Returns

tuple(float64[3N,3N], float64[3N,3N], float64[3N,3N], float64[3N,3N])

Some information

Notes

A super Matrix \(M=[[M1, M2], [M3, M4]]\) is invertible, if both the diagonal blocks, \(M_1\) and \(M_4\) are invertible The inverse of a (2x2) supermatrix can be calculated by [21], [22]:

(27)\[\begin{split}\begin{align*} & T_1 = (M_1 - M_2 M_4^{-1} M_3)^{-1} \\ & T_2 = -M_1^{-1} M_2 (M_4-M_3 M_1^{-1} M_2)^{-1} \\ & T_3 = -M_4^{-1} M_3 (M_1-M_2 M_4^{-1} M_3)^{-1} \\ & T_4 = (M_4 - M_3 M_1^{-1} M_2)^{-1} \\ \end{align*}\end{split}\]

pyrid.molecules.hydro_util.transform_reverse_supermatrix(zeta_tt_Tr, zeta_tr_Tr, zeta_rt_Tr, zeta_rr_Tr, pos, radii)

Transforms a (3N,3N) matrix into a (N,N,3,3) matrix.

Parameters

zeta_tt_Trfloat64[3N,3N]

Translational friction matrix.

zeta_tr_Trfloat64[3N,3N]

Translation-rotation coupling.

zeta_rt_Trfloat64[3N,3N]

Rotation-translation coupling.

zeta_rr_Trfloat64[3N,3N]

Rotational friction matrix.

Returns

tuple(float[N,N,3,3], float[N,N,3,3], float[N,N,3,3], float[N,N,3,3])

zeta_tt, zeta_tr, zeta_rt, zeta_rr

pyrid.molecules.hydro_util.transform_supermatrix(mu_tt, mu_tr, mu_rt, mu_rr, pos, radii)

Transforms a (N,N,3,3) matrix into a (3N,3N) matrix.

Parameters

mu_ttfloat64[N,N,3,3]

Translational mobility tensor.

mu_trfloat64[N,N,3,3]

Translation-rotation coupling.

mu_rtfloat64[N,N,3,3]

Rotation-translation coupling.

mu_rrfloat64[N,N,3,3]

Rotational mobility tensor.

Returns

tuple(float[3N,3N], float[3N,3N], float[3N,3N], float[3N,3N])

mu_tt_Tf, mu_tr_Tf, mu_rt_Tf, mu_rr_Tf