(255ap) Improving Computational Cost of Monte Carlo Models in RAFT Polymerization Processes

Reversible Addition-Fragmentation Transfer (RAFT) polymerization has become one of the most important types of Controlled Radical Polymerization (CRP) due to its potential application for the synthesis of new materials using a wide variety of monomers and mild reaction conditions.

Even though a substantial number of scientific publications have been directed to understanding the mechanism of the RAFT polymerization, there are still some questions about underlying details of the basic steps of the RAFT mechanism that remain unanswered. For example, no decisive explanation has been reported for the slowdown of the propagation rate of some RAFT systems when the concentration of the Chain Transfer Agent (CTA) is increased. Three main theories attempt to explain this rate retardation effect: the Slow Fragmentation (SF) theory, the Intermediate Radical Termination (IRT) theory and the Intermediate Radical Termination with Oligomers (IRTO) theory.

In this scenario, mathematical modeling becomes a valuable tool for helping to elucidate details of the reaction mechanisms. The time evolution of homogeneous chemical reactions is commonly calculated by solving a set of coupled ordinary differential equations. In this method, called the Deterministic Method, the species concentrations are represented by continuous functions of time and the predicted results are satisfactory in most cases. However, the formulation and solution of deterministic models usually involve some mathematical complexity and these models lack certain desired features that are important in polymerization processes. For instance, they do not take into consideration the stochastic nature of polymerization reactions. Besides, they are able to provide limited information about the microstructure of the polymer chains, when compared to other simulation methods.

The other important approach for modeling polymer processes is the stochastic simulation. This approach is mainly represented by the Monte Carlo (MC) technique, and specifically by the Gillespie¹ algorithm This algorithm, also called the Stochastic Simulation Algorithm (SSA), considers that the molecules are individually represented, is relatively simple to implement and offers extremely detailed information about the polymer microstructure and chain topology. This information can provide a valuable insight about the mechanistic details of the RAFT polymerization.

There are several different formulations of the SSA. Two of the most important are the so called Direct Method (DM) and the t-leaping method. The DM is characterized by simulating a single reaction in each iteration step. For this reason it is considered an exact simulation algorithm. However, it has the disadvantage of being slow and requiring large computational resources for some reaction systems. On the other hand, the t-leaping method divides the time history of the system into subintervals and calculates the number of times each reaction takes place in each subinterval. Hence, the method leaps along time from one subinterval to the next. The t-leaping method is one of the most used strategies in order to decrease the simulation time. The cost of this time reduction is a lower accuracy, compared to the DM. In this work we analyze implementations of the DM.

Regardless of the method used, the accuracy of these stochastic simulations depends significantly on the size of the sample of the reaction medium considered. The size of this sample, which can be defined by the number of molecules it contains, directly affects the computation times. The large values required for accurate and reproducible results can sometimes lead to prohibitive simulation times. Thus, the efficiency of the algorithm implementation becomes of critical importance, especially for the simulation of complex reaction systems.

In a previous work² a Monte Carlo model of a RAFT polymerization system was developed. It was implemented in Julia, a modern programming language designed to achieve high-performance in numerical and scientific computing. This model was used to obtain detailed information about the reaction system and the polymer molecular structure in RAFT polymerization processes. It provided an efficient method for predicting average properties and molecular weight distributions of all polymeric species present in the reaction system, including the bivariate MWD of the intermediate two-arm adduct. The model produced accurate and reproducible results of the MWDs in approximately 5 minutes of computing time on a regular desktop computer equipped with an Intel® Coreâ?¢ i5-3330 Processor (running at 3.20 GHz) and 8 Gb of RAM.

In the present work, different implementations of the DM algorithm are analyzed and compared for the simulation of RAFT reactions in order to improve the performance of our previous model. The algorithm is based on the work by Gillespie with several modifications to improve its time requirements. In one of the steps of this algorithm, all possible chemical reactions that can take place at a given instance during the simulation are listed, its reaction rates and occurrence probabilities calculated, and, based on these probabilities, one of the reactions is randomly chosen to be simulated next. This reaction selection step was found to be critical in the performance of the algorithm. Some strategies for dealing efficiently with this step have been proposed in the literature, such as the First Reaction Method³ (FRM), the Next Reaction Method⁴ (NRM), the Optimized Reaction Method⁵ (ORM) and the Sorting Direct Method⁶ (SDM). Most of these models were proposed for the simulation of biological processes, but they can easily be extended for polymerization reactions or other system of chemically reacting molecules.

In this work, we propose a modification of the ORM strategy to improve the efficiency of the reaction selection step of the DM algorithm. This strategy is based on the fact that in polymerization processes there are some specific reactions steps that occur considerably more frequent than others. The ORM algorithm takes advantage of this fact by placing the most frequent reactions at the beginning of the reaction search order. Since the MC simulation involves millions of iterations in which the reaction selection is repeated, a significant increase in the performance can be achieved. Unfortunately, in RAFT polymerization the sorting of chemical steps from most frequent to less frequent ones does not stay constant during the reaction time. This transient shift in reaction frequencies opens the possibility to further improve the performance of the ORM algorithm, by re-sorting the chemical reactions during the simulation and using specific data structures for reducing memory size and gaining speed performance.

The proposed model updates the sorted reaction array at different intervals based on the extent of reaction. This easy to compute quantity provides a reliable means to identify the shifts in the reaction frequency order. It does not depend on the reaction system individually and can be generalized for different mechanisms. Thanks to this modification and other minor optimizations the new model considerably improves the performance of our previous work.

(1) Gillespie, D. T. J. Phys. Chem. 1977, 81 (25), 2340â??2361.

(2) Pintos, E.; Sarmoria, C.; Brandolin, A.; Asteasuain, M. AIChE Annu. Meet. Proc. 2015.

(3) Gillespie, D. T. J. Comput. Phys. 1976, 22 (4), 403â??434.

(4) Gibson, M. A.; Bruck, J. J. Phys. Chem. A 2000, 104 (9), 1876â??1889.

(5) Cao, Y.; Li, H.; Petzold, L. J. Chem. Phys. 2004, 121 (9), 4059.

(6) McCollum, J. M.; Peterson, G. D.; Cox, C. D.; Simpson, M. L.; Samatova, N. F. Comput. Biol. Chem. 2006, 30 (1), 39â??49.