Next Article in Journal
Environmental Assessment of Trace Metals in San Simon Bay Sediments (NW Iberian Peninsula)
Previous Article in Journal
Application of Evolved Gas Analysis Technique for Speciation of Minor Minerals in Clays
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

KIMERA: A Kinetic Montecarlo Code for Mineral Dissolution

1
TECNALIA, Basque Research and Technology Alliance (BRTA), Parque Tecnológico de Bizkaia, Astondo Bidea, Edificio 700, 48160 Derio, Spain
2
Centro de Física de Materiales CSIC-UPV/EHU, Paseo Manuel de Lardizabal, 5, 20018 San Sebastian, Spain
3
Donostia International Physics Center, Paseo Manuel Lardizabal 3, 20018 San Sebastian, Spain
4
Department of Condensed Matter Physics, University of the Basque Country UPV/EHU, Barrio Sarriena s/n, 48940 Leioa, Spain
*
Author to whom correspondence should be addressed.
Minerals 2020, 10(9), 825; https://doi.org/10.3390/min10090825
Submission received: 3 August 2020 / Revised: 10 September 2020 / Accepted: 14 September 2020 / Published: 18 September 2020

Abstract

:
KIMERA is a scientific tool for the study of mineral dissolution. It implements a reversible Kinetic Monte Carlo (KMC) method to study the time evolution of a dissolving system, obtaining the dissolution rate and information about the atomic scale dissolution mechanisms. KIMERA allows to define the dissolution process in multiple ways, using a wide diversity of event types to mimic the dissolution reactions, and define the mineral structure in great detail, including topographic defects, dislocations, and point defects. Therefore, KIMERA ensures to perform numerous studies with great versatility. In addition, it offers a good performance thanks to its parallelization and efficient algorithms within the KMC method. In this manuscript, we present the code features and show some examples of its capabilities. KIMERA is controllable via user commands, it is written in object-oriented C++, and it is distributed as open-source software.

1. Introduction

Mineral dissolution has been widely studied due to its impact in important phenomena like soil formation [1,2], water and petroleum reservoir stability [3,4] or carbon sequestration [5,6,7,8], among others. The studies performed in the last decades have highlighted the importance of the nanoscale mechanisms in the dissolution process [9,10,11,12,13,14,15,16]. To verify and complement these experimental results, atomic scale computational methods like Density Functional Theory (DFT), molecular dynamics, molecular mechanics, Monte Carlo (MC) and Kinetic Monte Carlo (KMC) have been used by an increasing number of authors [17,18,19,20,21,22,23,24,25]. Of special interest is the KMC method which allows temporal scale studies comparable to experiments. The effect on dissolution of dislocations [26,27], grain sizes [25] or surface roughness [21], the inherent topographies associated to the dissolution mechanisms [22,23,26,27], the experimentally observed pulsating frequency at the nanoscale [28], or more recently the dissolution rate dependence with Δ G [29] are some of the milestones achieved by KMC simulations. Unfortunately, while many DFT, molecular dynamics and molecular mechanics programs are available both commercially and with free license [30,31,32,33], this is not the case with KMC. KMC simulations present the great advantage of being applicable in a multitude of fields, but also the disadvantage of being too specific to be programmed in a flexible and general package. SPPARKS [34], MONTECOFFE [35] and KMOS [36] are some KMC tools that were created to conduct studies in the field of materials science. Grain growth in annealed system [37], Mesoscale evolution in electron beam welding [38], adsorption of methane in zeolites [35] or oxidation of CO over Pt nanoparticles [35] are some examples of the studies performed with these codes. Nevertheless, this codes are not specific to the study of mineral dissolution, and performing dissolution simulations could be difficult or even not possible in certain systems. With that aim, we have developed KIMERA.
KIMERA, the name of which is an acronym for Kinetic Monte Carlo for Mineral Dissolution, seeks for being a helpful and efficient code which allows to any user with a basic knowledge on geochemistry to define and simulate the dissolution of a multitude of systems and minerals.

2. The Reversible Kinetic Monte Carlo Model for Mineral Dissolution

The KMC is an stochastic method based on the Markov property of independent processes happening within a system [39]. If the considered processes are chemical reactions, as it is the case of mineral dissolution, the Transition State Theory (TST) defines the reaction rate using only two parameters; the energy barrier of the reaction E B (J mol−1), and the frequency with which the system attemps to overcome that energy barrier f (s−1). Both parameters can be obtained from ab initio simulations. The rate r for a reaction from the initial state o to the final state p follows an Arrhenius equation:
r o p = f · exp E B k B · T
where k B is the Boltzmann constant ( 1.38 · 10 23 J K−1) and T the temperature (K).
The reactions, or events as commonly known in KMC method, happen sequentially in time with a probability proportional to the total rate r tot of all the reactions and a random number generation. This gives to the overall reaction the characteristic randomness of stochastic processes. The time increment to complete an event is also chosen randomly within a Poisson distribution of the total event probability:
Δ t = 1 r tot ln 1 u
where u is a random number uniformly distributed u ( 0 , 1 ] and r tot is the total rate of every single event:
r tot = o p o r o p o
p o accounts that each state o can access to different number of states. When an event happens, the system changes and all the possible transitions from the new state are redefined. Usually the transition affects only locally to a small part of the whole system, which can save great computational time [35]. This process is repeated until reaching a desired time or number of events (see Figure 1b).
The implementation on KIMERA is based on the reversible Kinetic Model that we presented recently [29]. In this model, we take into account microscopic reversibility defining for every possible dissolution reaction two events, one of dissolution and one of precipitation, with their respective rates r D and r P :
r D = f D · exp E D ( n ) k B · T
r P = f P · exp E P ( n ) Δ G * k B · T
where f i are the fundamental frequencies of each event, E D (n) and E P (n) the event energy barrier, and Δ G * the local free energy difference between the solid and the solution. Note that the precipitation expression includes the free energy difference between the solid and the solution. At far from equilibrium conditions, the Δ G * is large, and the precipitation term becomes negligible. Close to equilibrium Δ G * tends to zero, and the interplay between precipitation and dissolution reactions allows to mimic the dissolution rate decrease as the system reaches equilibrium.

3. The KIMERA Code

KIMERA aims to be broadly used either in big computer clusters or in personal computers. Therefore, we payed special attention to the implementation of user-friendly commands, a good performance and portability. It is written in the standard C++11, which ensures its portability. The program recognizes as input data a wide list of commands that can be found online https://[email protected]/mgp9999/kimera-publico.git. Thanks to these commands, the user can define the simulated system and obtain output files for visualization and analysis as well as restarting files. The program is based on the N-fold-way algorithm, which ensures a good efficiency [40]. Moreover, to overcome the typical problem of KMC simulations of waste of computer power in fast and repeating events [41], we have implemented a Poisson approximation to handle opposite effects of dissolution and precipitation which shows no reduction of performance or biased results [29]. If faster simulations were needed, KIMERA has been parallelized using openMP library [42], which entails an easy way to divide the workload of some loops of the program between the number of cores available.

3.1. Operation

The workflow of the program can be differentiated in three parts (see Figure 1a):

3.1.1. System Definition

The user defines the essential parameters of the simulation in an input file. The order in which commands are given is relevant, as some instructions can overwrite previous ones, totally, or partially. In these cases KIMERA always takes into account the last command. Some important steps are:
  • The mineral structure. The program can either read it from a standard ‘.xyz’ file, or can be defined by commands, or a combination of both. The ‘.xyz’ file [43] is easily obtained by tools such as VESTA [44] from downloadable ‘.cif’ files in mineralogical databases [45]. In principle KIMERA is thought to construct mineral surfaces replicating a small unit cell. Nevertheless, it is also possible to define a complex system within a ‘.xyz’ file and threat it as the whole system box. Coarse grained systems can be also simulated.
  • The system dimensions. The program replicates the unit cell in the three spatial directions. Studies of different planes are possible by unit cell transformations with external programs such as VESTA [46]. KIMERA can apply periodic boundary conditions (PBC) in the three spatial directions.
  • System shape. The program has commands to create different crystalline shapes of the system into complex surfaces or particles. For the moment the available geometries are cubic, spherical, parallelepiped, ellipsoidal, tick planes, or a combination of them.
  • Topographic defects can be defined in the system, such as insoluble regions, dislocations, impurities and vacancies. There are two ways of defining impurities; it is possible to define them in the unit cell indicating their occupancy, or introducing them ex post once the system has been defined.
  • Event definition. The KMC algorithm simulates the time evolution of a system as a set of possible events. These events take place at a rate that follows an Arrhenius equation (Equation (1)). A recent study demonstrates that the net dissolution of a mineral can be characterized using decoupled reactions of dissolution and precipitation [29]. Hence, we use that KMC scheme, so the fundamental frequency of the Arrhenius equation, f, splits into f D or f P , and E B into E D or E P depending if considering a dissolution or a precipitation event respectively. The energy barrier is characteristic of each chemical reaction and its neighbourhood, and can be obtained from the bibliography and/or ab initio calculations [22,23]. Supposing n neighbours of an atom, KIMERA can set E D and E P as a linear (Equation (6)) or a specific (Equation (7)) function of each neighbour j [23,47] (see Figure 2):
    E D k = E d · n E P k = E p · n
    E D k = f ( j ) E P k = g ( j )
    Note that Equation (6) is a specific case of Equation (7). Moreover, since the contribution to the energy barrier can be determined for several types of neighbours, k represents each set of contributors with the same characteristics and E D k and E P k its contribution for dissolution and precipitation energy barrier respectively.
    E D = k = 1 m E D k E P = k = 1 m E P k
    With these two ways of defining the energy barrier, two different approaches can be considered to describe the dissolution events:
    1
    A bond by bond description: Each linking bond breaks sequentially so that when an atom has no bonds left, it is released from the mineral.
    2
    A site by site description: All bonds reactions are unified in only one event, and each site dissolves with joint probability.
    As an additional element, KIMERA supports conditional event definition. Furthermore, it is possible to define the events based on ghost positions in the unit cell without physical meaning and to make a differentiation between atoms of the same type, for example it is possible to split the atoms of silica into Si1, Si2, etc. in the unit cell and then define events for each sub-element.
  • Target time (s). Predicting the time scale beforehand in a complex system can be tough. There are two options for the simulation to finish. The simplest option is to indicate the number of simulation steps, that is, the number of events to accomplish. The other option is to specify the target time (s) until the simulation is going to run. The user can request the program to do an estimation of it by considering the initial possible events. Given s initial sorted groups of rates corresponding to atom removals with different coordination r 1 < r 2 < < r s , the program approximates the total time for the system to dissolve as if all atoms N a had the same rate value; the previous to the middle one.
    t approx = N a r s 2 1
    This approximation arises from considering as limiting step the removal of the atom which leads to a kink atom, always with half of the mineral coordination [29,48].
  • Optional parameters related with the output files. As we will see, output files contain information of the system time evolution like snapshots for visualization or the quantity of dissolved atoms.

3.1.2. System Preparation

Once KIMERA has read the input commands, taking into account the event definition and the PBC, it elaborates a list with all the neighbours for each atom. This step consumes a relatively large computational time, and KIMERA can use the output file from previous simulations as starting point to improve the performance. If restarting files are used, some previous commands have necessarily to be the same, but some commands can be changed to modify the simulation conditions. Most common changes are the values of Δ G or the energy barriers.
From the neighbour data, the program elaborates a list with the reactive initial surface. Nevertheless, the surface can be modified ad-hoc by the user via commands.

3.1.3. Simulation Process

A key step of the KMC algorithm is to find a random event from the list of possible events looking at its rate [39] (see Figure 1b). This step is the most time consuming. KIMERA has two possible ways to do the search. By default, the binary search method [49] is done. Briefly, it consists of successively dividing the search range in two and discarding the one without the wanted value (see Figure 3b). The other method is the common lineal search, which compares one by one all the values of the list. The relation between the computing time and the number of elements to compute N is named computational complexity O. The complexity of binary search is O ( log ( N ) ) . The complexity of lineal search is higher, O ( N ) , but it can be parallelized. Later we will discus the performance of both search methods in Section 3.2 and Figure 3a.
During the simulation, the output files generated by KIMERA are:
  • Initial KIMERA file of the system in its own format (‘.initialkimerabox’). It is designed to save time in calculating neighbourhood, linked and affected atoms. A later simulation which reads this file will not need to do the preparation step.
  • Final KIMERA file of the system (‘.finalkimerabox’). When the simulation has finished, or has encountered an error, the system is printed in KIMERA format.
  • File with system snapshots (‘.box’) in LAMMPS format [33] for visualization. As this file can contain a lot of data, it may be better to handle the surface file unless for checking reasons.
  • File with surface snapshots (‘.surface’) in LAMMPS format. Instead of the whole system, only the atoms on the surface are printed in this file.
  • Data file with the time evolution of the following parameters (‘.data’): The total number of atoms dissolved of each type, its fraction, the surface dispersion, the gyradius (in no PBC systems) as well as all their derivatives.
  • Coordination file with the mean coordination to each type of atom along the dissolution process (‘.meandiscoord’). This data is key to calculate correctly Δ G value as explained below.
  • Layer atom files with the amount of atoms in each layer and each spatial direction (‘.alayer’, ‘.blayer’, ‘.clayer’). For example, the ‘.clayer’ file contains the total number of atoms of the cells in plane ab, layer by layer in c direction.

3.2. Parallelization Level

The parallelization level of a program is defined as the maximum speedup that a program can have. The speedup of a program from parallelization is limited by how much of the program can be parallelized [50]. For example, if 90% of the program can be parallelized, the theoretical maximum speedup using parallel computing would be 10 times no matter how many processors are used.
We have used the openMP library [42] to parallelize our code. As we have seen, the program presents three different parts: The system definition, the system preparation by KIMERA, and the simulation itself. The last two are the more time consuming and thus, they define the total simulation parallelization level. In Figure 3 the performance of the two first examples studied afterwards in Section 5 is plotted.
The preparation phase of the Kossel crystal example is very quick and hardly influences the total simulation time. Therefore it is a good example to get the parallelization level of the simulation phase. A decrease of 5 % by increasing the number of cores from 1 to 8 is obtained when using binary search. With lineal search the decrease is higher, 16 % . Such difference is due to the former search method cannot be parallelized. While lineal search seems to be more efficient, the roles are expected to be swapped in simulations with bigger systems and low number of cores.
On the other hand, the quartz grain example needs a long system preparation time, which has been tracked separately. While the total time presents a reduction of 14 % with 8 cores, the system preparation phase shows a good parallelization level with a reduction of 61 % . Therefore, best strategy to reduce the simulation time in a study with the same system is to use several cores to print only the initial Kimera file, to later used it in a set of subsequent simulations with only one processor and the default binary search.

4. Gibbs Free Energy Difference, Δ G

The Gibbs free energy difference Δ G is the driving force of a hydrolysis reaction in a mineral and it is related to the concentration of the chemical species of the mineral in the solvent. It is possible to obtain the dependence of the dissolution rate with Gibbs free energy in KMC simulations. The origin of this dependence has been demonstrated to reside at atomic scale and can be modeled by the interplay between dissolution and precipitation reactions, i.e., the microscopic reversibility [29]. The Gibbs free energy difference Δ G is related to the ion activity product, Π , of the dissolved material in water, divided by the solubility product K s :
Δ G = R · T · ln Π / K s = R · T · ln β
R is the ideal gas constant ( 8.3145 J mol−1 K−1) and T the temperature (K). The term Π / K s is more commonly known as saturation index β and sets the distance to equilibrium, where no dissolution happens.
Due to concentration gradients, and accumulation of anions and cations in the Stern layer [51], there is a difference between macroscopic concentrations in solution and that close the surface. Therefore, we can define a “local” or microscopic free energy difference Δ G * :
Δ G * = R · T · ln β *
In principle, the dissolution will be driven by this microscopic free energy difference. The local and macroscopic Δ G can then be related by:
Δ G = Δ G * R · T · ln β * β
Numerically, we cannot establish a direct relationship between both quantities, because we do not know the microscopic saturation index β * . However, we can relate them using the thermodynamic equilibrium condition and the dissolution and precipitation reaction rates. The thermodynamic equilibrium condition is define as the state when Δ G = 0 . In that state, the net reaction rate must be also equal to zero, and therefore r D = r P . Using these conditions, together with Equations (4) and (5), we can obtain the following relationship for a monocomponent mineral:
Δ G = Δ G * + E D j E P j ln f D f P
In Equation (13) we assume that the total dissolution and precipitation rates correspond to those of the reaction step that limits the dissolution close to equilibrium, usually a kink site. In addition, if a multicomponent mineral is considered, the macroscopic Δ G is given by coupling the concentration of each constituent element [48]. The difference between the macroscopic and local free energies is implicitly taken into account in the E D and E P values. Supposing a mineral with l atoms types, m sets of contributor, and n k neighbours for each contributor set, the Gibbs free energy difference, Δ G , is related with local Δ G * by considering that the net dissolution rate of the kink atoms of each type ( n ¯ k for all k) is 0 when Δ G = 0 ( k B T units):
Δ G = i = 1 l χ i · N ¯ i · Δ G * + k = 1 m E D k ( n ¯ k ) E P k ( n ¯ k ) ln f D f P
where χ i is the fraction of atoms of type i, N ¯ i is its average amount of broken bonds to dissolve, n ¯ k is the average number of neighbours of each contributor set to the atom type i in a bond breakage or formation, and E D ( n ¯ k ) and E P ( n ¯ k ) are their respective energy barrier values. Both χ i and n ¯ k can be obtained from the output data of KIMERA, from the ‘.data’ and ‘.meandiscoord’ files respectively. Since these values can change with the considered Δ G * value, the relation between Δ G * and Δ G may not be constant. Nevertheless, in practice the deviation is not high and can be considered as constant by calculating them from simulation at far from equilibrium conditions, when Δ G * and no precipitation events take place. The user must identify the value of N ¯ i by recognising the number of broken bonds during a dissolution process. This differentiation arises for example if we want to group several bond breaking events in only one event of different f and E B , such as a coarse grain simulation. In the examples of the following section the Δ G calculus is explicitly highlighted.

5. Examples of KIMERA Capabilities

5.1. Model A-B Kossel Crystals

In this section the dissolution with Δ G of two different configurations of a Kossel crystal consisting of two elements is described. A Kossel crystal, or Terrace Ledge Kink system (TLK), is a simple mineral structure consisting of a cubic lattice with six first neighbours [52]. Despite the simplicity of this system, it ensures enough topographic details so as to reproduce the mechanisms attributed to the dissolution process.
The systems in this model have two elements, A and B. In one of them, named ‘configuration 1’, the atoms are distributed alternately along x and y axes (not in z) within a box of 120 × 120 × 30 atoms (see Figure 4a). In the other one, named ‘configuration 2’ groups of 2 × 2 × 2 atoms of the same element are distributed alternately along x, y and z axes (see Figure 4b). In both cases one dislocation of 2 × 2 atoms wide and 20 atoms depth is set in the center of the system. The dissolution study is done in the {001} plane using PBC.
Both cases are similarly defined. The simulation parameters for the ‘configuration 2’ are first described, and then the minor changes for the other case are indicated.
  • The system dimensions are indicated; 60 × 60 × 15 unit cells.
  • The unit cell parameters. We used a = b = c = 5 Å and α = β = γ = 90 . Inside the cell, we define the 8 positions of the atoms in the unit cell that later on is repeated along the system. The positions are: (0,0,0), (2.5,0,0), (0,0,2.5), (2.5,0,2.5), (0,2.5,0), (2.5,2.5,0), (0,2.5,2.5) and (2.5,2.5,2.5) Å. All the positions are initially define as A atoms, and we later will redefine half of them as B type. Note that although the unit cell has 5.0 Å the distance between atoms is 2.5 Å, which is a typical distance reported for minerals [54].
    In ‘configuration 1’, the positions are the same, but half of them are of type B. Specifically, atoms in (0,0,0), (0,0,2.5), (2.5,2.5,0) and (2.5,2.5,2.5) are B type. Since the alternating disposition of the atoms is already taken into account with this definition, no additional commands to modify the system are needed.
  • We set periodic boundary conditions along x and y axes.
  • Physical parameters for the simulation. The target time t target = 4.0 · 10 2 s and the local Δ G * = 100 k B T units, which ensures far from equilibrium conditions. In ‘configuration 1’ the time scale is different due to its faster dissolution and t target = 4.0 · 10 1 s.
  • Event definition. We have chosen for this example an energy barrier for A-A atoms of E D A - A = 12 k B T units, for B-B E D B - B = 4 k B T units, which are respectively the higher and lower limit value for most minerals [27]. For the AB interaction the energy barrier is obtained from the Lorentz–Berthelot rule [55], E D A - B = E D A - A · E D B - B = 6.92 k B T units. The precipitation energy barrier for all the cases is E P A - B = E P A - A = E P B - B = 1 k B T units.
    For the frequency f D = f P = 4 · 10 13 s−1 s which lies in the range of values for atomic vibration in a mineral ( 10 10 10 13 s−1 at 300 K) [56].
    Lastly, KIMERA requires the number of neighbours that a bulk atom has to later define the initial reactive surface. For both for A and B atoms, a bulk atom has 3 A type neighbours and 3 B type neighbours. In ‘configuration 1’ the event definition is similar, but the number which defines a bulk atom changes. A bulk A atom has 2 A and 4 B neighbours. A bulk B atom has 4 A and 2 B neighbours.
  • Topographic defects. We define the last plane z = 0 as insoluble and we include one partial dislocation in the center with one third of the system depth. Since there are atoms within the dislocation that have a lower coordination than a bulk atom, the program recognised them as initial reactive surface. Therefore, we remove such condition because it is physically meaningless.
With those instructions KIMERA can start the simulations. Once the simulations are finished, we can relate Δ G * and Δ G with Equation (14). The ‘.meandiscoord’ file reports the average bond breakage of each atom type during all the simulation. The final steady values are collected in Table 1. Besides, the ‘.data’ files report the same dissolution rate for A and B atoms in both cases, thus χ A = χ B = 0.5 . The relation between Δ G * and Δ G is calculated from these values and Equation (14), and it is shown in Table 1. By changing the Δ G * value, the dissolution rate versus Δ G is obtained in the Figure 4f,g from the slope of the number of atoms removed versus time normalized to the surface area ( ( 60 · 5.0 · 10 10 ) 2 = 9 · 10 16 m2) and moles. Note that the rate may not be constant during the system dissolution. In ‘configuration 2’ the dissolution rate is constant with time since the dislocation has no effect. In contrast, in ‘configuration 1’ the dislocation opening provokes a progressive increase of the dissolution rate until dislocation coalescence. Thus we have taken the dissolution rate when it reaches the steady value after coalescence.
In ‘configuration 1’ the dislocation is not opened until enough energy is reached at Δ G = 9 kcal mol−1. The shape of the pit at this point is square (Figure 4c) but turns into a circular one (Figure 4d) with a small energy deviation. Besides, once the Δ G = 16.5 kcal mol−1 is reached, the spontaneous pit opening or mechanism III [29] is produced.
Regarding ‘configuration 2’, at close to equilibrium conditions only the very first B exposed atoms on the surface dissolve until only A atoms are exposed. The dissolution rate in this first zone is close to 0. Once the Δ G is negative enough at 27 kcal mol−1 the release of A atoms starts and the dissolution rate increases drastically. It presents a random dissolution pattern over all the surface without any influence of the dislocation (see Figure 4e). The onset is only due to the beginning of the dissolution of the limiting A atoms. Note that, as expected, with the same E D and E P energies the dissolution rate is much lower in ‘configuration 2’ since groups of A atoms restrict it. Moreover the energy needed to activate its dissolution is much higher; the onset is at much lower Δ G .
This type of theoretical examples on a Kossel crystal are very useful to study a general behaviour valid as a first approach for most minerals. However, as detailed below, more specific studies are possible in KIMERA.

5.2. Quartz Model I: Dissolution of an Ellipsoidal Grain

KIMERA is able to simulate the dissolution of real minerals with full atomistic structure resolution. As example, we simulate the dissolution of an ellipsoid grain of quartz with Δ G using the SCS-L1 model described by Kurganskaya and Luttge [23]. Briefly, the SCS-L1 model proposed in [23] considers that the energy barrier for a silicon atom dissolution depends linearly with the amount of first surrounding silicon atoms n 1 and, to a lesser extend, with the second surrounding silicon atoms n 2 and hydroxils groups n 3 bonded to the surface.
E B = n 1 · E B Si - Si 1 + ( n 2 + n 3 ) · E B Si - Si 2
Note that n 2 + n 3 = 12 in quartz. The specific energies are those proposed by [23]. Herein, we describe the instructions to run the simulation:
  • System dimensions. A box in which we will define the ellipsoid is created with 50 × 40 × 30 unit cells.
  • The unit cell parameters. For α -quartz a = b = 5.01 , c = 5.47 , α = β = 90 and γ = 120 . Inside the cell, we load a ‘.xyz’ file containing the positions, which has been converted from a ‘.cif’ file downloaded from The American Mineralogist crystal structure database [45]. Oxygen atoms can be removed for performance purposes since they are not explicitly taken into account for the quartz dissolution reaction in this case. The dissolution of a SiO2 is considered in a single step with a joint probability (Equation (6)). This can be interpreted as a coarse grain of a SiO2 unit in each Si site.
  • Physical parameters. The target time t target = 2.0 · 10 20 s and the local Δ G * = 100 k B T units.
  • Topographic defects. An ellipsoid with radius in the three axes, r x = 65 Å, r y = 85 Å and r z = 75 Å is defined as the simulation system. A dislocation along the x axis is placed in the middle.
  • Event definition. The energy barrier with first neighboring silicon atom is E D Si - Si 1 = 28 k B T units and with second E D Si - Si 2 = 4 k B T . Precipitation energies of E P Si - Si 1 = 10 k B T and E P Si - Si 2 = 1 k B T are used. All the four first silicon neighbours are at 3.09832 Å. If an atom is surrounded by the four first neighbours, it is considered to be in bulk. 12 second silicon neighbours are at 5.01 Å. Finally, the fundamental frequencies values are f D = f P = 1.0 · 10 12 s−1 [57].
The relation between Δ G * and Δ G is calculated using Equation (14) looking at the ‘.meandiscoord’ file to get the average breakage of atoms (see Table 2).
Figure 5 reports the dissolution rate versus Δ G and the time evolution of the grain, which is similar for all Δ G values. The grain dissolves maintaining an ellipsoidal shape with an irregular surface until its complete dissolution. The dissolution rate changes constantly as the exposed surface decreases. Therefore, we report the values when half of the forming atoms remains. The surface value is taken as a geometrical approximation in this point (see Figure 5b), A 1 / 2 = 4 · 10 16 m2.
The results show little influence of the dislocation in the dissolution rate. The local coordination at the dislocation does not decrease with respect to other grain regions, so it is not a preferential spot for dissolution. In this case the dissolution rate decreased gradually while getting close to equilibrium following a typical TST curve [56].

5.3. Quartz Model II: Dissolution of a Wulff Shape Particle

The previous quartz dissolution model is not able to reproduced jointly the experimental dissolution rate and activation energy [58]. If the dissolution parameters are fitted to reproduce the macroscopic activation energy, the dissolution rate is overestimated by several orders of magnitude, and vice versa. In this section we describe an example based on an alternative model proposed in [58]. In that work, in contrast to the previous site-by-site dissolution model, a bond-by-bond scheme is used, which is able to obtain both quantities along with the experimental topographies. The system is a wulff shape particle, the geometry with minimum surface energy of a mineral. For quartz it is a prism with hexagonal section and pyramidal tips [59] (see Figure 6a).
Instead of considering a coarse grain approximation where the quartz dissolves site by site, each oxygen bond breakage is explicitly considered. The energy barrier to break or reform the bond depends on its connectivity. We consider the connectivity by evaluating the number of unreacted oxygen atoms around the oxygen atom of interest. All of them are in a cutoff distance of 2.58 ± 0.01 Å. If one hydrolysis reaction occurs in one of this neighbouring linking oxygen atoms, the bond breaks reducing the connectivity and, consequently, the energy barrier for the hydrolysis of the oxygen atom of interest. This way it can be defined the energy barrier for a hydrolysis as a function of the surrounding oxygen atoms in bridging sites (see Table 3).
  • System dimensions. A box in which the wulff shape fit in is created with 16 × 16 × 47 units cells
  • Same unit cell parameters as the previous example. a = b = 5.01 , c = 5.47 , α = β = 90 and γ = 120 . The ‘.xyz’ file is also called, but this time the oxygen atoms do play an important role and they cannot be removed.
  • This time instead of target time, we define a target step S target = 62,700 steps, which is approximately the total amount of silicon and oxygen atoms forming the particle.
  • Topographic defects. The wulff shape is sculpted from the system by defining planes in which the atoms are removed. The equations of these planes are taken from GEODEBRA3D tool [60] which was used to visualise the system beforehand. Besides, two dislocations are defined and inner atoms removed from the initial surface. One dislocation is placed transversally in the center of the {100} plane, and another one perpendicular to the previous and longitudinally to the wulff shape
  • Event definition. The E D ( n ) and E P ( n ) for the linking oxygen bond is directly related with the number of oxygen atoms within 2.585 Å. 6 surrounding oxygen atoms indicates that the considered one is in a bulk position and therefore it is not reactive. Besides, the silicon atom must be automatically released if all of its four surrounding oxygen atoms have reacted. Finally, the fundamental frequency values f D = f P = 2.45 · 10 10 s−1 are indicated [61]. The Δ G * = 200 k B T value sets very far from equilibrium conditions.
The next step is to determine the relation between Δ G * and Δ G . From the ‘.meandiscoord’ file, most oxygen atoms react when three oxygen atoms are around ( n ¯ = 3 ). That means that the reference position to determine the Δ G value is a Q1-Q4 with E D = 24.1 and E P = 20.3 k B T units at 473 K (see Table 3). The relation between Δ G * and Δ G is Δ G = Δ G * + 3.8 .
It is important to highlight that in this process the number of broken bonds is N ¯ i = 1 , not like the previous model where the broken bonds corresponded to the number of surrounding silicon atoms.
Results of dissolution rate with Δ G is shown in Figure 6. They follow a typical TST curve. For all the simulations with Δ G , initial wulff shape is distorted into an elongated grain which reduces its size until the complete extinction. Same as the previous model, dislocations do not represent a preferential spot for dissolution. The surface area and the dissolution rate values are calculated at the time when the grain has lost half of its atoms (see Figure 6b). The estimated surface area is A 1 / 2 = 3.3 · 10 16 m2.

6. Conclusions

In this paper we have introduced KIMERA, an open source KMC code to study mineral dissolution. The code offers portability thanks to its implementation in C++11 language and a good performance thanks to efficient algorithms and its parallelization feature. KIMERA is a versatile code thanks to the multiple ways to define reaction events, and the possibility of performing studies with specific atomistic structures, coarse grained systems, periodic systems or particles.
To illustrate it, three different models in four different systems are studied. First, a lineal dependence of the energy barrier of dissolution with the number of first neighbours is considered to study two different Kossel crystal systems. Second, a lineal dependence of the energy barrier with the number of first and second neighbours is considered to study the dissolution of an ellipsoidal grain of quartz. Finally, a specific dependence of the energy barrier with surroundings to represent the sequential breakage of bridging oxygen bonds is used to study the quartz wulff shape dissolution.
Future enhancements of KIMERA will include: (1) An improvement and widening of the event definition, to even consider a differentiation of Δ G * , f D and f P with positions, (2) an extension of the code to consider growth and (3) the development of a tool, webpage or IDE to create and display the system of study.
KIMERA is available as open-source software under the GNU General Public License. Thus, KIMERA can be used free of charge, everyone can contribute to the software, extend it to his own needs and share newly developed plug-ins with other users. The C++ source code of KIMERA, all input file examples, as well as the ‘.xyz’ needed for the quartz examples can be downloaded from the bitbucket repository https://[email protected]/mgp9999/kimera-publico.git.

Author Contributions

Conceptualization, P.M., J.J.G., J.S.D. and H.M.; software, P.M.; validation, P.M., J.J.G., J.S.D. and H.M.; writing–original draft preparation, P.M.; writing–review and editing, P.M., J.J.G., J.S.D. and H.M. All authors have read and agreed to the published version of the manuscript.

Funding

J.S.D. acknowledges the funding from the Spanish Ministry of Economy, Industry and Competitiveness (project Ref-201860I057) and the Spanish Ministry of Science, Innovation and Universities (project Ref RTI2018-098554-B-I00). P. Martin acknowledges support from the PhD scholarship Tecnalia Research & Innovation’s grant.

Acknowledgments

All the simulations have been carried out at the high performance computing service of Basque Country i2basque and the UPV/EHU cluster. The authors thank for technical and human support provided by SGIker - UPV/EHU.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Van Breemen, N.; Buurman, P. Soil Formation; Springer: Berlin/Heidelberg, Germany, 2002. [Google Scholar]
  2. Kalbitz, K.; Solinger, S.; Park, J.H.; Michalzik, B.; Matzner, E. Controls on the dynamics of dissolved organic matter in soils: A review. Soil Sci. 2000, 165, 277–304. [Google Scholar] [CrossRef]
  3. Krouse, H.R.; Viau, C.A.; Eliuk, L.S.; Ueda, A.; Halas, S. Chemical and isotopic evidence of thermochemical sulphate reduction by light hydrocarbon gases in deep carbonate reservoirs. Nature 1988, 333, 415. [Google Scholar] [CrossRef]
  4. Canals, M.; Meunier, J.D. A model for porosity reduction in quartzite reservoirs by quartz cementation. Geochim. Cosmochim. Acta 1995, 59, 699–709. [Google Scholar] [CrossRef]
  5. Lal, R. Carbon sequestration. Philos. Trans. R. Soc. Lond. B Biol. Sci. 2008, 363, 815–830. [Google Scholar] [CrossRef] [PubMed]
  6. Lal, R. Soil Carbon Sequestration Impacts on Global Climate Change and Food Security. Science 2004, 304, 1623–1627. [Google Scholar] [CrossRef] [Green Version]
  7. Blain, S.; Quéguiner, B.; Armand, L.; Belviso, S.; Bombled, B.; Bopp, L.; Bowie, A.; Brunet, C.; Brussaard, C.; Carlotti, F.; et al. Effect of natural iron fertilization on carbon sequestration in the Southern Ocean. Nature 2007, 446, 1070. [Google Scholar] [CrossRef]
  8. Giammar, D.E.; Bruant, R.G., Jr.; Peters, C.A. Forsterite dissolution and magnesite precipitation at conditions relevant for deep saline aquifer storage and sequestration of carbon dioxide. Chem. Geol. 2005, 217, 257–276. [Google Scholar] [CrossRef]
  9. Juilland, P.; Gallucci, E. Morpho-topological investigation of the mechanisms and kinetic regimes of alite dissolution. Cem. Concr. Res. 2015, 76, 180–191. [Google Scholar] [CrossRef]
  10. Juilland, P.; Gallucci, E.; Flatt, R.; Scrivener, K. Dissolution theory applied to the induction period in alite hydration. Cem. Concr. Res. 2010, 40, 831–844. [Google Scholar] [CrossRef]
  11. Lasaga, A.C.; Blum, A.E. Surface chemistry, etch pits and mineral-water reactions. Geochim. Cosmochim. Acta 1986, 50, 2363–2379. [Google Scholar] [CrossRef]
  12. Brand, A.S.; Bullard, J.W. Dissolution kinetics of cubic tricalcium aluminate measured by digital holographic microscopy. Langmuir 2017, 33, 9645–9656. [Google Scholar] [CrossRef] [PubMed]
  13. Brand, A.S.; Feng, P.; Bullard, J.W. Calcite dissolution rate spectra measured by in situ digital holographic microscopy. Geochim. Cosmochim. Acta 2017, 213, 317–329. [Google Scholar] [CrossRef] [PubMed]
  14. Feng, P.; Brand, A.S.; Chen, L.; Bullard, J.W. In situ nanoscale observations of gypsum dissolution by digital holographic microscopy. Chem. Geol. 2017, 460, 25–36. [Google Scholar] [CrossRef] [PubMed]
  15. Cama, J.; Ganor, J.; Ayora, C.; Lasaga, C.A. Smectite dissolution kinetics at 80 °C and pH 8.8. Geochim. Cosmochim. Acta 2000, 64, 2701–2717. [Google Scholar] [CrossRef]
  16. Dove, P.M.; Han, N.; De Yoreo, J.J. Mechanisms of classical crystal growth theory explain quartz and silicate dissolution behavior. Proc. Natl. Acad. Sci. USA 2005, 102, 15357–15362. [Google Scholar] [CrossRef] [Green Version]
  17. Shvab, I.; Brochard, L.; Manzano, H.; Masoero, E. Precipitation mechanisms of mesoporous nanoparticle aggregates: Off-lattice, coarse-grained, kinetic simulations. Cryst. Growth Des. 2017, 17, 1316–1327. [Google Scholar] [CrossRef] [Green Version]
  18. Wang, Q.; Manzano, H.; Guo, Y.; Lopez-Arbeloa, I.; Shen, X. Hydration mechanism of reactive and passive dicalcium silicate polymorphs from molecular simulations. J. Phys. Chem. C 2015, 119, 19869–19875. [Google Scholar] [CrossRef]
  19. Wang, Q.; Manzano, H.; López-Arbeloa, I.; Shen, X. Water adsorption on the β-dicalcium silicate surface from DFT simulations. Minerals 2018, 8, 386. [Google Scholar] [CrossRef] [Green Version]
  20. Manzano, H.; Gartzia-Rivero, L.; Bañuelos, J.; López-Arbeloa, I.N. Ultraviolet–visible dual absorption by single BODIPY dye confined in LTL zeolite nanochannels. J. Phys. Chem. C 2013, 117, 13331–13336. [Google Scholar] [CrossRef]
  21. de Assis, T.A.; Reis, F.D.A. Dissolution of minerals with rough surfaces. Geochim. Cosmochim. Acta 2018, 228, 27–41. [Google Scholar] [CrossRef]
  22. Kurganskaya, I.; Luttge, A. Kinetic Monte Carlo approach to study carbonate dissolution. J. Phys. Chem. C 2016, 120, 6482–6492. [Google Scholar] [CrossRef]
  23. Kurganskaya, I.; Luttge, A. Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization. J. Phys. Chem. C 2013, 117, 24894–24906. [Google Scholar] [CrossRef]
  24. Rohlfs, R.D.; Fischer, C.; Kurganskaya, I.; Luttge, A. Crystal dissolution kinetics studied by a combination of Monte Carlo and Voronoi methods. Minerals 2018, 8, 133. [Google Scholar] [CrossRef] [Green Version]
  25. Briese, L.; Arvidson, R.S.; Luttge, A. The effect of crystal size variation on the rate of dissolution—A kinetic Monte Carlo study. Geochim. Cosmochim. Acta 2017, 212, 167–175. [Google Scholar] [CrossRef]
  26. Lasaga, A.C.; Luttge, A. Variation of crystal dissolution rate based on a dissolution stepwave model. Science 2001, 291, 2400–2404. [Google Scholar] [CrossRef]
  27. Meakin, P.; Rosso, K.M. Simple kinetic Monte Carlo models for dissolution pitting induced by crystal defects. J. Chem. Phys. 2008, 129. [Google Scholar] [CrossRef]
  28. Fischer, C.; Luttge, A. Pulsating dissolution of crystalline matter. Proc. Natl. Acad. Sci. USA 2018, 115, 897–902. [Google Scholar] [CrossRef] [Green Version]
  29. Martin, P.; Manzano, H.; Dolado, J.S. Mechanisms and Dynamics of Mineral Dissolution: A New Kinetic Monte Carlo Model. Adv. Theory Simul. 2019, 2, 1900114. [Google Scholar] [CrossRef] [Green Version]
  30. Artacho, E.; Anglada, E.; Diéguez, O.; Gale, J.D.; García, A.; Junquera, J.; Martin, R.M.; Ordejón, P.; Pruneda, J.M.; Sánchez-Portal, D.; et al. The SIESTA method; developments and applicability. J. Phys. Condens. Matter 2008, 20, 064208. [Google Scholar] [CrossRef]
  31. Gale, J.D. GULP: A computer program for the symmetry-adapted simulation of solids. J. Chem. Soc. Faraday Trans. 1997, 93, 629–637. [Google Scholar] [CrossRef]
  32. Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Nardelli, M.B.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys. Condens. Matter 2017, 29, 465901. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics; Technical Report; Sandia National Labs.: Albuquerque, NM, USA, 1993.
  34. Plimpton, S.; Thompson, A.; Slepoy, A. Stochastic Parallel PARticle Kinetic Simulator; Technical Report; Sandia National Laboratories: Albuquerque, NM, USA, 2008.
  35. Jørgensen, M.; Grönbeck, H. MonteCoffee: A programmable kinetic Monte Carlo framework. J. Chem. Phys. 2018, 149, 114101. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Hoffmann, M.J.; Matera, S.; Reuter, K. Kmos: A lattice kinetic Monte Carlo framework. Comput. Phys. Commun. 2014, 185, 2138–2150. [Google Scholar] [CrossRef] [Green Version]
  37. Holm, E.A.; Hoffmann, T.D.; Rollett, A.D.; Roberts, C.G. Particle-Assisted Abnormal Grain Growth; IOP Conference Series: Materials Science and Engineering; IOP Publishing: Risø, Denmark, 2015; Volume 89, p. 012005. [Google Scholar]
  38. Rodgers, T.M.; Madison, J.D.; Tikare, V.; Maguire, M.C. Predicting mesoscale microstructural evolution in electron beam welding. JOM 2016, 68, 1419–1426. [Google Scholar] [CrossRef] [Green Version]
  39. Jansen, A.P.J. An Introduction to Kinetic Monte Carlo Simulations of Surface Reactions; Springer: Berlin/Heidelberg, Germany, 2012; Volume 856. [Google Scholar]
  40. Bortz, A.B.; Kalos, M.H.; Lebowitz, J.L. A new algorithm for Monte Carlo simulation of Ising spin systems. J. Comput. Phys. 1975, 17, 10–18. [Google Scholar] [CrossRef]
  41. Dybeck, E.C.; Plaisance, C.P.; Neurock, M. Generalized temporal acceleration scheme for kinetic monte carlo simulations of surface catalytic processes by scaling the rates of fast reactions. J. Chem. Theory Comput. 2017, 13, 1525–1538. [Google Scholar] [CrossRef]
  42. Tian, X.; Bik, A.; Girkar, M.; Grey, P.; Saito, H.; Su, E. Intel® OpenMP C++/Fortran Compiler for Hyper-Threading Technology: Implementation and Performance. Intel Technol. J. 2002, 6, 36–46. [Google Scholar]
  43. Unofficial XYZ File Format Specification. Available online: http://en.wikipedia.org/wiki/XYZ_file_format (accessed on 1 October 2019).
  44. Momma, K.; Izumi, F. VESTA: A three-dimensional visualization system for electronic and structural analysis. J. Appl. Crystallogr. 2008, 41, 653–658. [Google Scholar] [CrossRef]
  45. Downs, R.T.; Hall-Wallace, M. The American Mineralogist crystal structure database. Am. Mineral. 2003, 88, 247–250. [Google Scholar]
  46. Unit Cell Rotation with VESTA. Available online: https://ma.issp.u-tokyo.ac.jp/en/app-post/1788 (accessed on 15 March 2020).
  47. Kohli, C.; Ives, M. Computer simulation of crystal dissolution morphology. J. Cryst. Growth 1972, 16, 123–130. [Google Scholar] [CrossRef]
  48. Lasaga, A.C.; Lüttge, A. Mineralogical approaches to fundamental crystal dissolution kinetics. Am. Mineral. 2004, 89, 527–540. [Google Scholar] [CrossRef]
  49. Williams, L.F., Jr. A modification to the half-interval search (binary search) method. In Proceedings of the 14th Annual Southeast Regional Conference, Birmingham, AL, USA, 22–24 April 1976; pp. 95–101. [Google Scholar]
  50. Hill, M.D.; Marty, M.R. Amdahl’s law in the multicore era. Computer 2008, 41, 33–38. [Google Scholar] [CrossRef] [Green Version]
  51. Dukhin, S.S.; Kretzschmar, G.; Miller, R. Dynamics of Adsorption at Liquid Interfaces: Theory, Experiment, Application; Elsevier: Amsterdam, The Netherlands, 1995. [Google Scholar]
  52. Oura, K.; Katayama, M.; Saranin, A.; Lifshits, V.; Zotov, A. Surface Science; Springer: Berlin/Heidelberg, Germany, 2003; Volume 602, pp. 229–233. [Google Scholar] [CrossRef]
  53. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. Model. Simul. Mater. Sci. Eng. 2010, 18, 015012. [Google Scholar] [CrossRef]
  54. Gibbs, G.V.; Downs, R.T.; Cox, D.F.; Ross, N.L.; Prewitt, C.T.; Rosso, K.M.; Lippmann, T.; Kirfel, A. Bonded interactions and the crystal chemistry of minerals: A review. Z. Krist. Cryst. Mater. 2008, 223, 1–40. [Google Scholar] [CrossRef]
  55. Berthelot, D. Sur le mélange des gaz. Compt. Rendus 1898, 126, 1703–1855. [Google Scholar]
  56. Lasaga, A.C. Kinetic Theory in the Earth Sciences; Princeton University Press: Princeton, NJ, USA, 2014. [Google Scholar]
  57. Pelmenschikov, A.; Leszczynski, J.; Pettersson, L.G. Mechanism of dissolution of neutral silica surfaces: Including effect of self-healing. J. Phys. Chem. A 2001, 105, 9528–9532. [Google Scholar] [CrossRef]
  58. Martin, P.; Gaitero, J.J.; Dolado, J.S.; Manzano, H. A Kinetic Monte Carlo Model for Quartz Dissolution. 2020; in preparation. [Google Scholar]
  59. Wendler, F.; Okamoto, A.; Blum, P. Phase-field modeling of epitaxial growth of polycrystalline quartz veins in hydrothermal experiments. Geofluids 2016, 16, 211–230. [Google Scholar] [CrossRef] [Green Version]
  60. Geodebra3d. Available online: https://www.geogebra.org/3d?lang=en (accessed on 15 March 2020).
  61. Casey, W.H.; Lasaga, A.C.; Gibbs, G. Mechanisms of silica dissolution as inferred from the kinetic isotope effect. Geochim. Cosmochim. Acta 1990, 54, 3369–3378. [Google Scholar] [CrossRef]
Figure 1. (a) KIMERA (Kinetic Monte Carlo for Mineral Dissolution) workflow. (b) KMC (Kinetic Monte Carlo) workflow.
Figure 1. (a) KIMERA (Kinetic Monte Carlo for Mineral Dissolution) workflow. (b) KMC (Kinetic Monte Carlo) workflow.
Minerals 10 00825 g001
Figure 2. Example of lineal or specific definition of energy barrier available in KIMERA. The orange sphere i represents the particle under consideration, the blue and red particles j represent particles bonded to the i of two different types. The sub-indexes of these j particles indicate their type k 1 and k 2 , and their number, 4 and 2 for k 1 and k 2 respectively.
Figure 2. Example of lineal or specific definition of energy barrier available in KIMERA. The orange sphere i represents the particle under consideration, the blue and red particles j represent particles bonded to the i of two different types. The sub-indexes of these j particles indicate their type k 1 and k 2 , and their number, 4 and 2 for k 1 and k 2 respectively.
Minerals 10 00825 g002
Figure 3. (a) Normalized simulation time of the examples with the number of cores. The A-B Kossel crystal example in blue and green (‘configuration 1’ explained below) presents with both search algorithms a simulation time with 1 core of 18885 s in our computer. The total simulation time for the quartz grain in red 2255 s with 1 core, and its preparation phase in orange 902 s. (b) Binary search working scheme. The target value is successively limited if it is not found.
Figure 3. (a) Normalized simulation time of the examples with the number of cores. The A-B Kossel crystal example in blue and green (‘configuration 1’ explained below) presents with both search algorithms a simulation time with 1 core of 18885 s in our computer. The total simulation time for the quartz grain in red 2255 s with 1 core, and its preparation phase in orange 902 s. (b) Binary search working scheme. The target value is successively limited if it is not found.
Minerals 10 00825 g003
Figure 4. A-B Kossel Crystal studies with Δ G . (a) Initial ‘configuration 1’ system. (b) Initial ‘configuration 2’ system. (c) Topography of ‘configuration 1’ at the beginning of the onset (red square). (d) Topography of ‘configuration 1’ at far from equilibrium conditions (blue point). (e) Topography of ‘configuration 2’ at far from equilibrium conditions (green point). (f) Dissolution rate versus Δ G for ‘configuration 1’. (g) Dissolution rate versus Δ G for ‘configuration 2’. The visual representation is done using OVITO program [53].
Figure 4. A-B Kossel Crystal studies with Δ G . (a) Initial ‘configuration 1’ system. (b) Initial ‘configuration 2’ system. (c) Topography of ‘configuration 1’ at the beginning of the onset (red square). (d) Topography of ‘configuration 1’ at far from equilibrium conditions (blue point). (e) Topography of ‘configuration 2’ at far from equilibrium conditions (green point). (f) Dissolution rate versus Δ G for ‘configuration 1’. (g) Dissolution rate versus Δ G for ‘configuration 2’. The visual representation is done using OVITO program [53].
Minerals 10 00825 g004
Figure 5. Quartz grain dissolution study with Δ G . (a) Initial grain. (b) Grain after half of the forming atoms have dissolved. (c) 10% of the atoms remain. (d) Dissolution rate versus Δ G curve. Blue point corresponds to the simulation of the topographies in (ac), though similar topographies are obtained in any other point. The visual representation is done using OVITO program.
Figure 5. Quartz grain dissolution study with Δ G . (a) Initial grain. (b) Grain after half of the forming atoms have dissolved. (c) 10% of the atoms remain. (d) Dissolution rate versus Δ G curve. Blue point corresponds to the simulation of the topographies in (ac), though similar topographies are obtained in any other point. The visual representation is done using OVITO program.
Minerals 10 00825 g005
Figure 6. Quartz wulff shape dissolution study with Δ G . (a) Initial state. (b) Wulff shape after half of the forming atoms have dissolved. (c) 25% of the atoms remain. (d) Dissolution rate versus Δ G curve. Blue point corresponds to the simulation of the topographies in (ac), though similar topographies are obtained in any other point. The visual representation is done using OVITO program.
Figure 6. Quartz wulff shape dissolution study with Δ G . (a) Initial state. (b) Wulff shape after half of the forming atoms have dissolved. (c) 25% of the atoms remain. (d) Dissolution rate versus Δ G curve. Blue point corresponds to the simulation of the topographies in (ac), though similar topographies are obtained in any other point. The visual representation is done using OVITO program.
Minerals 10 00825 g006
Table 1. Average bond breakage for each type of atom in the A-B Kossel crystal configurations and Δ G and Δ G * relation.
Table 1. Average bond breakage for each type of atom in the A-B Kossel crystal configurations and Δ G and Δ G * relation.
A-AA-BB-BB-A Δ G
Configuration 11.01.331.02.66 3 · ( Δ G * + 18.76 )
Configuration 21.50.4231.52.576 3 · ( Δ G * + 19.52 )
Table 2. Average breakage value for each bond for the quartz grain example and Δ G * and Δ G relation.
Table 2. Average breakage value for each bond for the quartz grain example and Δ G * and Δ G relation.
Si-Si-3.09832Si-Si-5.0100Si-Si-5.66774Si-Si-4.42416 Δ G
1.951.921.921.95 1.95 · ( Δ G * + 52.47 )
Table 3. Model energy barrier values with the bond surroundings in kJ mol−1.
Table 3. Model energy barrier values with the bond surroundings in kJ mol−1.
Bond E D E P Surrounding Oxygen Atoms
Q1-Q1--0
Q1-Q275601
Q1-Q385702
Q1-Q495803
Q2-Q285702
Q2-Q395803
Q2-Q41051054
Q3-Q31051054
Q3-Q41351355
Q4-Q4--6

Share and Cite

MDPI and ACS Style

Martin, P.; Gaitero, J.J.; Dolado, J.S.; Manzano, H. KIMERA: A Kinetic Montecarlo Code for Mineral Dissolution. Minerals 2020, 10, 825. https://doi.org/10.3390/min10090825

AMA Style

Martin P, Gaitero JJ, Dolado JS, Manzano H. KIMERA: A Kinetic Montecarlo Code for Mineral Dissolution. Minerals. 2020; 10(9):825. https://doi.org/10.3390/min10090825

Chicago/Turabian Style

Martin, Pablo, Juan J. Gaitero, Jorge S. Dolado, and Hegoi Manzano. 2020. "KIMERA: A Kinetic Montecarlo Code for Mineral Dissolution" Minerals 10, no. 9: 825. https://doi.org/10.3390/min10090825

APA Style

Martin, P., Gaitero, J. J., Dolado, J. S., & Manzano, H. (2020). KIMERA: A Kinetic Montecarlo Code for Mineral Dissolution. Minerals, 10(9), 825. https://doi.org/10.3390/min10090825

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop