=Paper= {{Paper |id=Vol-2267/170-174-paper-31 |storemode=property |title=Molecular dynamic simulation of water vapor interaction with various types of pores using hybrid computing structures |pdfUrl=https://ceur-ws.org/Vol-2267/170-174-paper-31.pdf |volume=Vol-2267 |authors=Vladimir V. Korenkov,Eduard.G. Nikonov,Maria. Popovičová }} ==Molecular dynamic simulation of water vapor interaction with various types of pores using hybrid computing structures== https://ceur-ws.org/Vol-2267/170-174-paper-31.pdf

Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018

MOLECULAR DYNAMIC SIMULATION OF WATER VAPOR
INTERACTION WITH VARIOUS TYPES OF PORES USING
HYBRID COMPUTING STRUCTURES
V.V. Korenkov1,3, a, E.G. Nikonov1, b, M. Popovičová2, с
1
Joint Institute for Nuclear Research, 141980 Dubna, Moscow Region, Russian Federation
2
University of Prešov, 080 01 Prešov, Slovakia
3 Plekhanov Russian University of Economics, 36 Stremyanny per., Moscow, 117997, Russia

E-mail: a korenkov@jinr.ru, b e.nikonov@jinr.ru, с maria.popovicova@unipo.sk

Theoretical and experimental investigations of water vapor interaction with porous materials are very
needful for various fields of science and technology. Mathematical modelling plays an important role
in these investigations. Conventional approaches to solve problems of mathematical research of the
processes of interaction of water vapor with individual pore are the so called micro and macro
approaches. The first approach is based on the use of diffusion equation for description of interaction
of water vapor with a pore. The second approach is based on various particle methods like Monte-
Carlo simulations, molecular dynamics (MD) etc. In MD usage of efficient calculation methods is
necessary because the degree of approximation for simulating system is largely determined by the
dimensionality of the system of equations being solved at every time step. Number of time steps is
also quite large. In this work, a study of efficiency of various implementations algorithms for MD
simulation of water vapor interaction with individual pore is carried out. It is also investigated
dependence of time required for simulations on different parameters, like number of particles in the
system, shape of pores, and so on. The results of parallel calculations are compared with the results
obtained by serial calculations.

Keywords: porous media, molecular dynamics, macroscopic diffusion model, parallel calculations

170
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018

1. Introduction
One of the most important problem in numerical simulation based on molecular dynamics or
Monte-Carlo approach of many particle systems is the need to use huge computing resources to obtain
realistic simulation results. A system of water vapor and a pore is an example of such many particle
systems. Theoretical and experimental investigations of water vapor interaction with porous materials
are very needful for various fields of science and technology. Not only studies of the interaction of
water vapor and porous material as a continuous medium, but also the study of the interaction of water
vapor with individual pore is very important in these researches. Mathematical modelling occupies an
important place in these investigations. Conventional approaches to solve problems of mathematical
research of the processes of interaction of water vapor with individual pore are the following.
The first approach is based on the use of diffusion equation for description of interaction of
water vapor with a pore. It is so called macro approach. The second approach is based on various
particle methods like, for example, molecular dynamics (MD). These methods essentially consider the
micro-structure of the investigated system consisting of water vapor and a pore. This second approach
can be called a micro approach.
At the macro level, the influence of the arrangement structure of individual pores on the
processes of water vapor interaction with porous material as a continuous medium is studied. At the
micro level, it is very interesting to investigate the dependence of the characteristics of the water vapor
interaction with porous media on the geometry and dimensions of the individual pore. Both
approaches require the most efficient calculation methods as far as possible with the current level of
development of computational technologies. Usage of efficient calculation methods is necessary
because the degree of approximation for simulating system is largely determined by the
dimensionality of the system of equations being solved at every time step. Number of time steps is
also quite large.
In this work, a study of efficiency of various implementations algorithms for MD simulation
of water vapor interaction with individual pore is carried out. A great disadvantage of MD is its
requirement of a relatively large computational effort and long time in simulations. These problems
can be drastically reduced by parallel calculations. In this work we investigate dependence of time
required for simulations on different parameters, like number of particles in the system, shape of
pores, and so on. The results of parallel calculations are compared with the results obtained by serial
calculations. Two-dimensional and three-dimensional models of the pore are used for comparative
analysis of parallel and serial calculations.

2. Molecular dynamics model
In classical molecular dynamics, the behavior of an individual particle is described by the
Newton equations of motion [1], which can be written in the following form
𝑑 2 𝑟⃗𝑖 (
𝑚𝑖 2 = 𝑓⃗𝑖 , 1)
𝑑𝑡
where 𝑖 − a particle number, 1 ≤ 𝑖 ≤ 𝑁, 𝑁 − the total number of particles, 𝑚𝑖 − particle
mass, 𝑟⃗𝑖 − coordinates of position, 𝑓⃗𝑖 − the resultant of all forces acting on the particle. This resultant
force has the following representation
𝜕𝑈(𝑟⃗1 , … , 𝑟⃗𝑁 ) (
𝑓⃗𝑖 = − + 𝑓⃗𝑖𝑒𝑥 ,
𝜕𝑟⃗𝑖 2)
where 𝑈 − the potential of particle interaction, 𝑓⃗𝑖 − a force caused by external fields. For a
𝑒𝑥

simulation of particle interaction, we use the Lennard-Jones potential [2] with 𝜎 = 3.17Å and
𝜀 = 6.74 ∙ 10−3 eV. It is the most used to describe the evolution of water in liquid and saturated vapor
phase. Equations of motion (1,2) were integrated by Velocity Verlet method [3]. Berendsen thermostat
[4] is used for temperature calibration and control. The coefficient of the velocity recalculation 𝜆(𝑡) at
every time step 𝑡 depends on the so called ''rise time'' of the thermostat 𝜏𝐵 which belongs to the
interval [0.1,2] ps. 𝜏𝐵 describes strength of the coupling of the system to a hypothetical heat bath. For
increasing 𝜏𝐵 , the coupling weakens, i.e. it takes longer to achieve given temperature 𝑇0 from current

171
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018

temperature 𝑇(𝑡). The Berendsen algorithm is simple to implement and it is very efficient for reaching
the desired temperature from far-from-equilibrium configurations. Initial concentrations were obtained
from the density of water vapor at the appropriate pressure and density at a given temperature using
known tabulated data. The pressure in the pore was controlled using the formula based on virial
equation [5].

1
𝑃 = (〈2𝐾〉 − ⟨∑ 𝑟𝑖𝑗 ∙ 𝑓(𝑟𝑖𝑗 )⟩ ).
3𝑉
𝑖<𝑗

Here 𝑉 is the pore volume, 〈2𝐾〉 is the doubled kinetic energy averaged over the ensemble,
𝑓(𝑟𝑖𝑗 ) is the force between particles 𝑖 and 𝑗 at a distance 𝑟𝑖𝑗.

3. Computational algorithm for molecular dynamic simulation
For molecular dynamic simulation we used the code written in CUDA C. The program does
not require a lot of memory. We only keep co-ordinates, speeds and forces for each particle. One of
the main problems of molecular dynamic simulation is a large number of particles and time steps.
Therefore, it is necessary to use parallel calculations. The code for our simulations was implemented
on heterogeneous computing cluster HybriLIT.
The code contains four functions that are paralleled, and which are performed on the GPU.
This is a function for calculating the forces (i.e., acceleration) for individual particles, which calculates
the interactions between all particles (F1). There are two functions to calculate new coordinates and
speeds for each particle. We need two functions to calculate them because we need forces acting on
particles at two different time moments (F2 and F3). Finally, we use the Berendsen thermostat in the
program that runs parallel at the GPU too (F4). Other calculations are performed on the host. General
scheme of the calculation algorithm for two- and three-dimensional molecular dynamic simulation is
shown in Figure 1.
In this paper we compare the temporal realization of these four functions F1-F4 on the GPU
and the CPU. The total time of parallel computing consists of two parts, that is the time needed
directly to calculate on the GPU (pure GPU time) and the time needed to complete these calculations
on the CPU because some algorithms performed on the GPU must be completed by the CPU. In this
work, total GPU time will indicate the sum of these two times.

4. 2D molecular dynamic simulation
For two-dimensional (2D) molecular dynamic simulation we consider 2D model of the pore in
the shape of a rectangle with dimensions 𝑙𝑥 = 1𝜇𝑚, 𝑙𝑦 = 1𝜇𝑚. The outer space in this micro-model
reflects is a rectangle right to the pore with 𝑘 ∙ 𝑙𝑥 and (2𝑘 + 1) ∙ 𝑙𝑦 dimensions. All sides of the outer
space satisfy to the periodic boundary conditions. The left pore side reflects the inner molecules due to
the boundary condition [6] but also provides the periodic boundary conditions for a part of outer
space.
There are 420 molecules of water vapor inside the pore which form saturated water vapor at
temperature 35℃ and pressure 5.62 kPa at the time t=0. The value of parameter k = 3 means that the
outer space volume for calculations is 21 times larger than the pore volume. There are 1764 molecules
of water vapor in outer space corresponding to 20% saturated water vapor. The integration step is
0.016 ps and evolution time is 65.3 ns.

172
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018

Figure 1. Parallelized algorithm for molecular dynamic simulation

5. 3D molecular dynamic simulation
In three-dimensional (3D) case [7] we made simulation for 3D model of the pore in the shape
of a prism of dimensions 𝑙𝑥 = 500𝑛𝑚, 𝑙𝑦 = 50𝑛𝑚, 𝑙𝑧 = 50𝑛𝑚. Five walls are isolated and so there is
no exchange of particles with outer space. The sixth wall is open and connected with the external
environment. The external environment is modelled by a prism which is 9 times bigger than the pore.

173
Proceedings of the VIII International Conference "Distributed Computing and Grid-technologies in Science and
Education" (GRID 2018), Dubna, Moscow region, Russia, September 10 - 14, 2018

The big prism satisfies periodic boundary conditions. This means that the particles which pass through
one wall return to the system through the opposite wall. The integration time step is 0.016 ps and
evolution time is 65.3 ns.
Consequently, we have considered the following input data for the drying process. There are
1000 molecules of water vapor inside the pore which forms saturated water vapor at temperature 25℃
and pressure 3.17 kPa. There are 1800 molecules of water vapor in the outer area space corresponding
to 20% saturated water vapor.

6. Conclusion
For both simulations, in 2D and 3D cases, the average time 𝑡1 required to calculate one step
was calculated. To compare the two simulations, we have converted this time to one particle 𝑡𝑝 . The
results are shown in Table 1. The time required to calculate one step for one particle for 2D and 3D
simulation on GPU in both cases (pure GPU time and total GPU time) especially pure GPU time is
2:3. CPU time does not keep this ratio.
Table 1. Comparison of computation times per step and for one particle for both simulations
2D 3D
Time
CPU pure GPU total GPU CPU pure GPU total GPU
𝑡1 112.732 4.532 4.668 184.349 8.649 8.786
𝑡𝑝 0.05162 0.00207 0.00214 0.06584 0.00309 0.00314

As our investigations showed for both cases of 2D and 3D simulation, when paralleling the
computations, there are some optimal value of number of threads in blocks 27 such that the
computation time becomes minimal then for other values of this number of threads. Parallelized
algorithm runs about 24 times faster than the sequential one in the 2D case and about 21 times faster in
the 3D case for one time step 𝑡1 and for average time step per particle 𝑡𝑝 . In addition, it should be
noted that, when parallelizing, the cost ratio of the computation time per particle for 2D and 3D
modeling is equal 2/3 with high precision.

References
[1] Gould H., Tobochnik J. and Christian W. An Introduction to Computer Simulation Methods. 3rd
edition. // Addison-Wesley, 2006. - 796 pp.
[2] Lennard-Jones J. E. On the Determination of Molecular Fields. // Proc. Roy. Soc. 1924. V. A106.
P. 463.
[3] Verlet L. Computer experiments on classical fluids. I. Thermodynamical properties of Lennard-
Jones molecules. // Phys. Rev. 1967. V. 159. P. 98.
[4] Berendsen H. J. C. et al. Molecular dynamics with coupling to an external bath. // J. Chem. Phys.
1984. V. 81. P. 3684.
[5] Frenkel D. and Smith B. Understanding molecular simulation: from algorithms to applications.
Second edition. // Academic Press, 2006. - 658 pp.
[6] Nikonov E.G., Pavluš M., Popovičová M. 2D microscopic and macroscopic simulation of water
and porous material interaction. 2017. Available at: arXiv:1709.05878 [physics.flu-dyn]
[7] Nikonov E.G., Pavluš M., Popovičová M. Molecular dynamic simulation of water vapor
interaction with blind pore of dead-end and saccate type. 2017. Available at: arXiv:1708.06216
[physics.flu-dyn]

174