Design of the Algorithm for Packaging of Water Molecules in a Fixed Volume

The question about the closest symmetrical packing of chemical substance species (molecules, ions, polymer chains, nanoparticles, etc.) is a subproblem of predicting the structure of matter, in particular, the structure of a crystal, information about which makes it possible to predict almost all of its properties. The design of mathematical models for the closest symmetrical packing is an important and a challenging task for the practical application of optimization theory in theoretical chemistry. In this paper, we report the development of the algorithm for water molecules’ symmetrical packing in a closed space of the rectangular parallelepiped. In our model, water molecules were represented as regular tetrahedrons, in which two vertices correspond to hydrogen atoms in the molecule (vertices of type A), and the other two correspond to electron pairs (vertices of type B). The connections between them are shown in Figure 1. It is also important to note that due to the chemical properties of water, bonds between molecules in our model can only be formed between hydrogen and an electron pair, i.e., only bonds of the form A∼B or B∼A will be considered. Thus, it is necessary that the algorithm considers the types of vertices when attaching molecules and does not allow connections of the form A∼A and B∼B. Note also that this type of calculation has been performed non-quantum mechanically, and no specific quantum effects (e.g., dealing with electron correlation) were taken into account in our model.

The relevance of the topic is supported by numerous papers on structure prediction [1]. In [2] K. K. Fisher studied a way to search for crystal structures using machine learning and modern quantum mechanical methods. N. Braun and G. Huttner developed an algorithm that shows correlations between molecular shape and the type of packing [3]. In [4] N. W. Thomas suggested an analytical approach to characterize organic molecular crystal structures where close packing is an important factor.

The research team from the State University of New York at Buffalo presented the evolutionary algorithm XtalOpt and its refinements [5]. A. R. Oganov and C. W. Glass developed a reliable methodology for predicting crystal structures based on evolution algorithms and have demonstrated the efficiency of this method for systems with up to 30–40 atoms per cell [6]. Then, A. O. Lyakhov, A. R. Oganov, H. T. Stokes, and C. Zhu Q. further refined the previous algorithm and achieved a reliable representation for systems of 200 atoms per cell [7]. It is worth noting that the algorithm presented in this work is capable of packing a larger number of atoms: up to 3000 in a rectangular parallelepiped.

E. V. Peresypkina and V. A. Blatov determined the crystal structures of 33,575 monosystem organic compounds within the first three coordination spheres, and demonstrated that generally molecular packings are constructed according to one of the close packings [8]. In [9], A. Gavezzotti and L. Lo Presti introduced a new view of crystal packing that results from joining static structural information from X-ray analysis with dynamic upgrades. The main advantage of this algorithm is that it allows us to quickly and efficiently estimate the number of molecules in a limited volume represented as a parallelogram. In the presented algorithm, a mathematical approximation of the water molecule is made, which makes it easy to obtain the structure of the arrangement of molecules in the parallelepiped. This approximation simplifies the calculations and reduces the running time of the algorithm, which further leads to the possibility of scaling the problem. At the same time, the possibility of adequate representation of the chemical structure of the substance is kept [10].

The designed algorithm has a practical value because it allows one to determine the number of water molecules, which can be placed in a rectangular parallelepiped with specified values of length, width and height. This algorithm can be easily used for the creation of initial solvent massive for molecular dynamics’ simulations and theoretical studies of the closest packing of molecules. Note that, in the future, we plan to adapt and extend this approach to both the sphere and the ellipsoid, which could be very relevant for fundamental colloid chemistry (problem of correct description of water drops). To the best of our knowledge, to date there are no MATLAB implementations for such algorithms in the literature. Hence, the constructed algorithm can be used as a new tool to estimate the densest packing in specified area. The algorithm consists of the following steps:

The cycle on step 8 of the algorithm uses a special procedure. It rotates the adding molecule to be connected so as to create a hydrogen–electron pair connection. This procedure determines the coordinates of the vertices of the tetrahedron corresponding to the added molecule, depending on the location of the previous molecule (subject to the types A∼B, B∼A communication requirements). Futhermore, this procedure takes into account the distance between molecules on the x, y, and z axes, which is set by the user, thus avoiding overlapping and adjusting the distance between molecules depending on the problem.

We implement the algorithm described above in MATLAB. The program calculates the number of molecules that could be stored in a given volume. Furthermore, it forms a table with the coordinates of the vertices of the tetrahedra in the Cartesian coordinate system centred at one of the vertices of the rectangular parallelepiped. This table can be used to interface the results with a visualization program, which required the coordinates of the hydrogen atoms to build atom-based system. The information about the vertices type (having hydrogen or not) can be easily added to this table, because the table stores information about tetrahedron vertices in an ordered way: of every four verticies of tetrahedron, the first two contain hydrogen and the other two do not. Moreover, the program has the ability to determine the interval between the molecules at any of the coordinates (by parameters $h_x,h_y,h_z$).

The MATLAB procedure is presented below. Firstly, we define the supporting functions.

1. Addition a tetrahedron to a graph

 

 

 

 

2. Determination of coordinates of the new tetrahedron along with the current vertex type and vertex coordinates with the largest abscissa

 

 

The body of the main program:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

We developed the symmetrical packing algorithm for water molecules and implemented it in the MATLAB environment. In addition, a series of numerical experiments were carried out. The advantage of the implemented algorithm is that, unlike the common problem of the dense packing of balls, it works with tetrahedra. In addition, algorithm takes into account chemical bonds and is capable of laying more than a thousand molecules in less than 6 s, while packing almost 4 thousand molecules in just 45 s (see Table 1). However, packing 8.5 thousand molecules will take more than 3 min. Thus, the algorithm’s running time increased almost 4.5 times after doubling the number of molecules. Therefore, for practical use, it is advisable to break down the entire volume so that the 4 thousand molecules are the largest number of molecules in each cell. Note also that this type of calculation has been performed non-quantum mechanically, and no specific quantum effects (e.g., deal with electron correlation) were taken into account in our model.

In order to illustrate the location of the water molecules obtained as a result of the algorithm, the program displays a graph with the image of tetrahedra (see Figure 2 and Figure 3). The asterisks mark the vertices of type A.