rigidbody_util

Short description

RBs

class pyrid.molecules.rigidbody_util.RBs(*args, **kwargs)[source]

A rigid body/bead model object is defined by a set of particles with fixed positions relative to the body center. In our framework we use rigid bodies as a course grained representation of molecules. A rigid bead model has a center of diffusion around which the rigid bead model rotates in response to torque and rotational diffusion. By default a surface molecule is bound to the mesh surface at the center of diffusion or shifted by an offset passed by the user. A rigid bodies topology is described by a list of particle indices. A RBs instance is initializeed without any molecules. Use the add_RB method to add molecules to the list.

Attributes

Dataarray_like: Numpy structured array containing all data that define a moelcules state. dtype = np.dtype([(‘next’, np.uint64),(‘name’, ‘U20’),(‘id’, np.int32),(‘type_id’, np.int32), (‘pos’, np.float64, (3,)), (‘dX’, np.float64, (3,)), (‘force’, np.float64, (3,)), (‘torque’, np.float64, (3,)), (‘topology’, np.int32, (20,)),(‘topology_N’, np.int32),(‘q’, np.float64, (4,)), (‘dq’, np.float64, (4,)),(‘B’, np.float64, (4,4)),(‘orientation_quat’, np.float64, (3,3)),(‘mu_tb’, np.float64, (3,3)),(‘mu_rb’, np.float64, (3,3)),(‘mu_tb_sqrt’, np.float64, (3,3)),(‘mu_rb_sqrt’, np.float64, (3,3)),(‘Dtrans’, np.float64),(‘Drot’, np.float64),(‘radius’, np.float64), (‘loc_id’, np.int32), (‘compartment’, np.int64), (‘triangle_id’, np.int32), (‘pos_last’, np.float64, (3,)), (‘Theta_t’, np.float64, (3,)), (‘Theta_r’, np.float64, (3,)), (‘posL’, np.float64, (3,)), (‘collision_type’, np.int64), (‘next_transition’, np.float64), ], align=True)
posarray: position of the center of mass

Methods

add_RB(name, compartment, System, Particles, vol_surf_id)	Adds a new rigid body molecule to the Data array.
next_um_reaction(System, i)	Calculates the time of the next unimolecular reaction.
set_triangle(tri_id, i)	Set the triangle id of a surface molecule.
add_particle(pi, i)	Adds a particle to the molecule topology list.
update_topology(System, Particles, i)	Updates the positions of all particles in the molecule.
place_particles(pos, types, radii, System, Particles, i)	Registeres the particle positions, radii and types of a molecule.
set_pos(Particles, System, pos, i)	Sets the position of moelcule i to pos.
update_force_torque(Particles, i)	Updates the total force and torque acting on molecule i.
update_B(i)	Updates matrix B which is used to update the moelcule’s rotation quaternion.
calc_orientation_quat(i)	Calculated the orientation matrix of the molecule from it’s’ rotation quaternion.
set_orientation_quat(Particles, q, System, i)	Rotates the molecule according to quaternion q.
update_orientation_quat(i)	Updates the molecules orientation matrix and the B-matrix.
update_dq(System, i)	Updates the rotational displacement dq.
update_dX(System, i)	Updates the translational displacement dX.
update_particle_pos(Particles, System, i)	Updates the position and orientation of volume molecule i (and all its particles).
update_particle_pos_2D(System, Particles, i)	Updates the position and orientation of surface molecule i (and all its particles).
rescale_pos(Particles, System, mu, i)	Rescales the position molecule i’s origin by a factor mu (needed for the Berendsen barostat).

__init__()[source]

add_RB(name, compartment, System, Particles, vol_surf_id)[source]

Adds a rigid bead molecule to the system.

Parameters

namestring: name of the molecule type
compartmentint64: Id of the compartment the moelcule sites in
Systemobject: Instance of System class
Particlesobject: Instance of Particles class

add_particle(pi, i)[source]

Adds a new particle to the rigid bead molecule topology.

Parameters

piint64: Particle index
iint64: Molecule index

Raises

IndexError(‘Maximum number of particles per molecule (20) reached! To change this, currently you have to manualy set the maximum size in rigidbody_util.py’): by default, currenztly, a maximum number of 20 particles per rigid bead molecule type is allowed. The user may incerase this limit by increasing the size of the field ‘topology’ in the structured numpy array dtype ‘item_t_RB’ (see rigi_body_util.py)

calc_orientation_quat(i)[source]

Calculates the orientation/rotation matrix of moelcule i from it’s rotational quaternion q.

Parameters

iint64: Molecule index

Notes

The rotation matrix corresponding to a unit rotation quaternion \(\boldsymbol{q} = [q_0, q_1, q_2, q_3]\) can be calculated by:

(1)\[\begin{split}\boldsymbol{A} = \begin{pmatrix} 2(q_0^2+q_1^2)-1 & 2(q_1 q_2-q_0 q_3) & 2(q_1 q_3-q_0 q_2) \\ 2(q_1 q_2-q_0 q_3) & 2(q_0^2+q_2^2)-1 & 2(q_2 q_3-q_0 q_1) \\ 2(q_1 q_3-q_0 q_2) & 2(q_2 q_3-q_0 q_1) & 2(q_0^2+q_3^2)-1 \\ \end{pmatrix}. \end{split}\]

(Note: There exist several equivalent ways to calculate the rotation matrix, which share the way the off-diagonal elements are calculated but slightly differ in the diagonal elements)

next_um_reaction(System, i)[source]

Draws the time point of the next unimolecular diffusion for moelcule i.

Parameters

Systemobject: Instance of System class
iint64: Molecule index

Notes

For unimolecular reactions we can draw the time point of the next reaction from a distribution using a variant of the Gillespie Stochastic Simulation Algorithm (SSA) Stiles and Bartol [11], Erban et al. [12]. For a single molecule having \(n\) possible transition reactions/reaction paths each having a reaction rate \(k_i\), let \(k_t = \sum_i^n k_i\) be the total reaction rate. Then, we can draw the time point of the next reaction from:

\[\tau = \frac{1}{k_t} \ln\Big[\frac{1}{r}\Big],\]

where \(r\) is a random number uninformly distributed in (0,1).

place_particles(pos, types, radii, System, Particles, i, h_membrane)[source]

Places all the particles that make up the topology of rigid bead molecule i given a list of partcile types, positions and radii.

Parameters

posfloat64[:,3]: Array of particle positions
typesarray_like: Array of strings of particle type names
radiifloat64[:]: Array of particle radii
Systemobject: Instance of System class
Particlesobject: Instance of Particles class
iint64: Molecule index
h_membranefloat64: Additive factor by which the particle positions are shifted along the z axis of the local coordinate frame. Thereby, the relative position to the center of diffusion is shifted. ‘h_membrane’ should only be used to change the position at which a surface molecule sits in the mesh surface.

rescale_pos(Particles, System, mu, i)[source]

Rescales the position of molecule i by a multiplicative factor mu. This method is only used for the Berendsen barostat.

Parameters

Particlesobject: Instance of Particles class
Systemobj: Instance of System class.
mufloat64: Multiplicative factor by which to scale the molecule position.
iint: Index of the molecule.

Returns

dtype: Some information

Raises

NotImplementedError (just an example): Brief explanation of why/when this exception is raised

set_orientation_quat(Particles, q, System, i)[source]

Sets the orientation of rigid bead molecule i given a rotation quaternion q.

Parameters

Particlesobject: Instance of Particles class
qfloat64[4]: Rotation quanternion
Systemobject: Instance of System class
iint64: Molecule index

set_pos(Particles, System, pos, i)[source]

Sets the position of rigid bead molecule i.

Parameters

Particlesobject: Instance of Particles class
Systemobject: Instance of System class
posfloat64[3]: Position of the molecule center
iint64: Molecule index

set_triangle(tri_id, i)[source]

Sets the triangle id for surface molecule i.

Parameters

triangle_idint64: Index of the triangle
iint64: Molecule index

update_B(i)[source]

Calculates matrix B used in the quaternion propagator (see RBs.dq()).

Parameters

iint64: Molecule index.

Notes

When propagating the cartesian coordinates of the rigid bead moelcule’s center of diffusion we need to account for the orientation of the moelcule. This is readily done by applying the rotation matrix A to the molecule’s mobility tensor (Note: to rotate a tensor we need to apply the rotation matrix twice: \(\mu^{lab} = A \cdot \mu^{local} \cdot A^{-1}\)). However, things are more complicated when deriving the propagator for the rotation quaternion. Here, the rotation matrix is partially replaced by matrix \(\boldsymbol{B}\) which accounts for the fact that we are using a four dimensional coordinate set of quaternions instead of three euler angles. In other terms, \(\boldsymbol{B}\) relates the angular velocity in the local frame to the quaternion velocity ( \(\frac{d\boldsymbol{q}}{dt} = \boldsymbol{B} \boldsymbol{\omega}^{local}\) ) [13]. The matrix \(\boldsymbol{B}\) is given by:

(2)\[\begin{split}\boldsymbol{B} = \frac{1}{2} \begin{pmatrix} q_0 & -q_1 & -q_2 & -q_3 \\ q_1 & q_0 & -q_3 & q_2 \\ q_2 & q_3 & q_0 & -q_1 \\ q_3 & -q_2 & q_1 & q_0 \\ \end{pmatrix},\end{split}\]

where \(\boldsymbol{q}\) is the current rotational quaternion of molecule i.

update_dX(System, i)[source]

Calculates the translational displacemnt of molecule i using the Brownian dynamics algorithm for anisotropic particles introduced by Ilie et al. [13].

Parameters

Systemobj: Instance of System class.
iint: Index of the molecule.

Notes

The translational displacement of an anisotropic particles due to Brownian motion and an external force \(\boldsymbol{F}\) can be expressed by [13]:

(3)\[\begin{split}\begin{align*} R_{\alpha}(t+\Delta t) - R_{\alpha}(t) = & A_{\alpha \gamma} \mu_{\gamma \delta}^{tb} A_{\beta \delta} F_{\beta} \Delta t \\ & + A_{\alpha \gamma}(\sqrt{\boldsymbol{\mu}^{tb}})_{\gamma \beta} \Theta_{\beta}^t \sqrt{2 k_B T \Delta t}, \end{align*}\end{split}\]

where \(\boldsymbol{A}\) is the rotation matrix and \(\boldsymbol{\mu}^{tb}\) the translational mobility tensor of the molecule in the body fixed frame (as opposed to the lab frame of reference). \(\boldsymbol{\Theta}^t\) is a normal distributed random vector describing the random movement of the molecule due to collisions with the fluid molecules. \(\boldsymbol{\mu}^{tb}\) should be a real, positive semidefinite matrix to have proper physical meaning [14] . In this case, there exists a unique square root of the matrix \(\sqrt{\boldsymbol{\mu}^{tb}}\), which can be found via diagonalization (see geometry_util.Transform:Util.sqrtM()). Latin indices run from 0 to 3 and Greek indices run from 1 to 3.

update_dq(System, i)[source]

Calculates the rotational displacemnt of molecule i using the Brownian dynamics algorithm for anisotropic particles introduced by Ilie et al. [13].

Parameters

Systemobj: Instance of System class.
iint: Index of the molecule.

Notes

The rotational displacement of an anisotropic particles due to Brownian motion and an external torque \(\boldsymbol{T}\) can be expressed by [13]:

(4)\[\begin{split}\begin{align*} q_{a}(t+\Delta t) - q_{a}(t) = & B_{a \alpha} \mu_{\alpha \beta}^{rb} A_{\gamma \beta} T_{\gamma} \Delta t \\ & + B_{a \alpha}(\sqrt{\boldsymbol{\mu}^{rb}})_{\alpha \beta} \Theta_{\beta}^q \sqrt{2 k_B T \Delta t} + \lambda_q q_a, \end{align*}\end{split}\]

where \(\boldsymbol{A}\) is the rotation matrix and \(\boldsymbol{\mu}^{rb}\) the rotational mobility tensor of the molecule in the body fixed frame (as opposed to the lab frame of reference). \(\boldsymbol{\Theta}^q\) is a normal distributed random vector describing the random rotational movement of the molecule due to collisions with the fluid molecules. \(\boldsymbol{\mu}^{rb}\) should be a real, positive semidefinite matrix to have proper physical meaning [14]. In this case, there exists a unique square root of the matrix \(\sqrt{\boldsymbol{\mu}^{rb}}\), which can be found via diagonalization (see geometry_util.Transform:Util.sqrtM()). Latin indices run from 0 to 3 and Greek indices run from 1 to 3.

The matrix \(\boldsymbol{B}\) is given by:

(5)\[\begin{split}\boldsymbol{B} = \frac{1}{2q^4} \begin{pmatrix} q_0 & -q_1 & -q_2 & -q_3 \\ q_1 & q_0 & -q_3 & q_2 \\ q_2 & q_3 & q_0 & -q_1 \\ q_3 & -q_2 & q_1 & q_0 \\ \end{pmatrix}. \end{split}\]

For infinitesimal time steps \(\boldsymbol{q}\) preserves its unit length. However, for finite time step \(\Delta t\) an additional constraint needs to be added to ensure that \(\boldsymbol{q}\) keeps its unit length. This is realized here by a constraint force in the last term, directed along \(\boldsymbol{q}\). By solving the Langrange multiplier \(\lambda_q\) from the condition \(q(t+\Delta t) = 1\) we get the strength of the force. Denoting \(\tilde{q}(t+ \Delta t)\) as the uncopnstrained quaternion, \(\lambda_q\) is obtained from solving the quadratic equation

(6)\[\lambda_q^2 + 2 \lambda_q \boldsymbol{q} \cdot \tilde{q}(t+ \Delta t)+ \tilde{q}^2(t+ \Delta t) = 1\]

Simple rescaling of the quaternion will change the sampled phase space distribution. However, in practice, the resulting error is marginal and since rescaling is fatser than solving the quadratic equation each time step, we currently do the rescaling. However, the more accurate method is mentioned here for compeleness.

update_force_torque(Particles, i)[source]

Updates the total force and torque acting on rigid bead molecule i by summing over all forces and torques of its particles.

Parameters

Particlesobject: Instance of Particles class
iint64: Molecule index

update_orientation_quat(i)[source]

Updates the rotation matrix of rigid bead molecule i, however, without updating the molecule and particle positions.

Parameters

iint64: Molecule index

update_particle_pos(Particles, System, i)[source]

Updates the position of volume molecule i.

Parameters

Particlesobject: Instance of Particles class
Systemobj: Instance of System class.
iint: Index of the molecule.