update_force

Short description

@author: Moritz F P Becker

pyrid.system.update_force.react_interact_test(comp_i, loc_type_i, mol_idx_j, mol_idx_i, RBs, System)

Tests if two particles are valid interaction or reaction partners.

Parameters

comp_iint64

Compartment index of molecule i

loc_type_iint64

Location type of moelcule i (volume : 0, surface : 1)

mol_idx_jint64

Index of moelcule j

mol_idx_iint64

Index of moelcule i

RBsobject

Instance of RBs class

Systemobject

Instance of System class

Notes

We need to consider several cases:

1. If both particles are in the same compartment

   1. and at least one is a volume molecule, interactions and reactions are always allowed.

   2. If both particles are surface molecules, interactions are always valid but reactions are only allowed if the angle between the molecules’ triangle normals is less than 90 degree.

2. If both particles are in different compartments

   1. and both particles are volume molecules, neither interactions nor reactions are allowed.

   2. If both particles are surface particles, interactions are allowed, however, not reactions except binding reactions!

   3. If one of the particles belongs to a surface molecule and the other to a volume molecule

      1. reactions are always allowed if one of the compartments is the simulation box (System), otherwise neither reactions, nor interactions are allowed.

We neglect reactions for surface molecules at an angle above 90 degree because the error introduced by using euclidian distances becomes large. In general, you should not use meshes with high local curvature (relative to the molecule size and reaction readius) and/or large angles between neighboring triangle for this reason. Also, triangles that are far apart in terms of the geodesic distance can come close to each other in terms of euclidian distance. In the latter case, the angle between the triangle is, however, in most cases larger than 90 degree. Thereby, neglecting reactions also prevents this type of erroneous reactions.

pyrid.system.update_force.update_force_append_reactions(Particles, System, RBs, HGrid)

Updates the pair interaction forces between the particle pairs and between particles and the compartemnt meshes (and the simulation box boundary in case of repulsive boundary conditions). Also adds the bimolecular and unimolecular particle reactions to the reactions list.

Parameters

Systemobj

Instance of System class.

Particlesobj

Instance of Particle class.

RBsobj

Instance of RBs (Rigid Bodies) class.

HGridarray_like

Hierarchical grid that divides the simulation box into cells. The cell size decreases with the level of the hierarchical grid. Particles are assigned to a level depending on their interaction cutoff radius.

dtype = np.dtype([(‘L’, np.int64, (N_Levels,)), (‘sh’, np.float64, (N_Levels,3)), (‘head’, np.int64, (np.sum(NCells)*2,)), (‘head_start’, np.int64, (N_Levels,)), (‘ls’, np.int64, (2*N+1,)), (‘NCells’, np.int64, (N_Levels,)), (‘cells_per_dim’, np.int64, (N_Levels,3)), (‘cutoff’, np.float64, (N_Levels,))], align=True)

Notes

Hierarchical Grid

The computationaly most expensive part in molecular dynamics simulations is usually the calculation of the pairwise interaction forces, beacuase to calculate these, we need to determine the distance between particles. When doing this in the most naiv way, i.e. for each particle i we iterate over all the other particles j in the system and calculate the distance, the computation time will increase quadratically with the number of partciles in the system (\(O(N^2)\)). Even if we take into account Newtons third law, this will decrease the number of computations (\(O(\frac{1}{2}N(N-1))\)), but the computational complexity still is quadratic in N. A straight forward method to significantly improve the situation is the linked cell list approach [3] (pp.195: 5.3.2 Cell structures and linked lists) where the simulation box is divided into \(n \times n \times n\) cells. The number of cells must be chosen such that the side length of the cells in each dimension \(s = L/n\), where \(L\) is the simulation box length, is greater than the maximum cutoff radius for the pairwise molecular interactions (and in our case also bimolecular reactions). This will decrease the computational complexity to \(O(14 N \rho s^3 )\), where \(\rho\) is the molecule density (assuming a mostly homogeneous distribution of molecules). Thereby, the computation time increase rather linear with \(N\) instead of quadratically (\(N^2\)).

However, one problem with this method is that it does not efficiently handle polydisperse particle size distributions. This becomes a problem when doing minimal coarse graining of proteins and other structures we find in the cell such as vesicles. As mentioned above, n should be chosen such that the cell length is greater than teh maximum cutoff radius. If we would like to simulate, e.g. proteins (\(r \approx 2-5 nm\)) in the presence of synaptic vesicles (\(r \approx 20-30 nm\)), the cells become way larger than necessary from the perspective of the small proteins.

One way to approach this problem would be, to choose a small cell size (e.g. based on the samllest cutoff radius) and just iterate not just over the nearest neighbour cells but over as many cells such that the cutoff radius of the larger proteins/SVs is covered. This approach has has a big disadvantages: Whereas for the smaller partciles we actually reduce the number of unnecessary distance calculations, for the larger particles we don’t. We can, however, still take advantage of Newton’s 3rd law. For this, we only do the distance calculation if the radius of particle i is larger than the radius of particle j. If the radii are equal, we only do the calculation if index i is smaller than index j.

A much better approach has been introduced by Ogarko and Luding [4] and makes use of a so called hierarchical grid. This approach is the one we use in PyRID. In the hierarchical grid approach, each particle is assigned to a different cell grid depending on its cutoff radius, i.e. the grid consists of different levels or hierarchies, each having a different cell size. This has the downside of taking up more memory, howewer, it drastically reduces the number of distance calculations we need to do for the larger particles and also takes advantage of Newtons third law, enabling polydisiperse system simulations with almost no performance loss. The algorithm for the distance checks works as follows:

1. Iterate over all particles

2. Assign each particle to a cell on their respective level in the hierarchical grid

3. Iterate over all particles once more

4. Do a distance check for the nearest neigbour cells on the level the current particle sites in. This is done using the classical linekd cell list algorithm.

5. Afterwards, do a cross-level search. For this, distance checks are only done on lower hierarchy levels, i.e. on levels with smaller partcile sizes than the current one. This way, we do not need to double check the same particle pair (3rd law). However, in this step, we will have to iterate over potentialy many empty cells (Note: this should, in principle, also be doable in a more efficient way).

Mesh collision response

The distance between particle and mesh triangles is calculated and based on the distance a repulsive force is determined. The repulsive force is normalized by the total number of triangle collisions. This approach is not exact but works sufficiently well. Particle vertex collisions are neglected if one or more particle face collisions have been detected.