GPU-Accelerated Large-Scale Molecular Simulation Toolkit
Station: Research

1. Introduction:

Coarse-grained (CG) models based on molecular dynamics (MD) from atomistic simulations are often utilized to describe polymers on the mesoscale. However, the main problem is that most available CG force fields do not treat bond breaking and bond formation, thus a correct description of the bonding process is still inaccessible. How to link polymer dynamics with polymerization kinetics on larger scale is still a challenge, which demands new computational simulation methods in which chemical reactions are taken into account. Notably, the “reaction” or “polymerization” in the following denotes more accurately a simulated chain growth process in CG level, but not a real chemical reaction process commonly used in quantum chemistry. To solve the problem of coupling MD-driven dynamics with “reactions” on the CG level, a model of treating both the particle dynamics and the reactive processes with force fields is developed by Hong Liu, You-Liang Zhu and co-workers. Generally, the CG reaction models can be categorized into two classes. For reproducing the chain propagation in polymerization, some works used many local reaction steps with the active growth centers connecting to free monomers within their capture radii. This class of technique could be identified as pure “distance-based” protocols, in which they let the active ends randomly search for the closest free monomer to connect to, so that the chains could grow. A typical free radical polymerization reaction is very fast, thus in a CG simulation the reaction turns out to be a pure probability issue of judging whether the reaction can take place within one time step or not. The simulation studies based on this idea can be identified as the second class, that is, “probability-based” protocols. Within the capture radius, the bond formation in the second class in each step is not a certain event, instead it is probability-based. The model here we proposed combines both protocols by providing two selectable reaction rules, i. e. reaction with the closest free monomer at a probability or the free monomer within the capture radius.

2. Model and methods:

where [M] is the free monomer concentration, [P*] is the concentration of growth centers in the system, and T is the reaction time interval.

3. Usage

3.1 Data format

The information of particles related to polymerization model could be read from XML file in following formats.

(1) The indication of initiator particles with int type, 1 for initiator.

$\langle$h_init num="4"$\rangle$






(2) The crosslinking number of particles with int type, 0 usally for reactable monomer.

$\langle$h_cris num="4"$\rangle$






3.2 General rules

(1) The polymerization models have three types including free-radical, step-growth, and exchange ones which are listed in following section. The polymerization model is free-radical type by default and could be implicitly changed to step-growth or exchange ones by setMaxCris() or setExchangePr() function, respectively. The polymerization model also could be explicitly called by FrpMode(), setSgapMode(), and setExchMode() functions for free-radical, step-growth, and exchange types, respectively.

(2) The particles with h_init=1 are active center for reaction. For free-radical and exchange polymerization models, the particles with h_cris=0 are reactive monomers, whareas particles with h_cris=1 are inert units for reaction. For step-growth polymerization model, the particles with h_cris < max_cris where max_cris with a default value of 1 could be changed by setMaxCris() function are reactive monomer. A reaction happens between active center and reactive monomer. Reaction will add 1 to the h_cris of active center and reactive monomer, respectively.

(3) The reaction between two active centers (both with h_init=1) is forbidden, unless function setInitInitReaction(True) is called. The reaction probability which could be set by setPr() for free-radical and step-growth polymerization model and setExchangePr() for exchange polymerization model consider the sequence of the parameters. The function setPr('A', 'B', 0.005) means that the particle A is active center, particle B is reactive monomer, and their reaction probability is 0.005, whereas function setPr('B', 'A', 0.005) means B is active center and particle A is reactive monomer. Similarly, the function setExchangePr ('A', 'B', 'C', 0.1) means particle A is reactive monomer, particle B is active center, and partice C is the particle connecting to particle B and will be replaced by particle A in a reaction.

3.3 Polymerization models

(1) Free-radical polymerization

Formula: AAAA* + A -> AAAAA*

Description: The polymerization starts from the particle A*, which represents an active center. The active center reacts with a reactive monomer A. After the reaction, the reactive center moves ahead to the new connected particle. Thereby, the chain grows with continous reactions with reactive monomers.

A script example:

reaction =galamost.Polymerization(all_info, neighbor_list, 2**(1.0/6.0), 16361)

reaction.setPr('A', 'A', 0.005)



Note: The functions reaction.setPr('A', 'A', 0.005) and reaction.setPeriod(50) indicate that there is a reacting probability of 0.005 every 50 time steps between active center A (initally with h_init=1) and reactive monomer A (with h_init=0 and h_cris=0). The fucntion setMaxCris('A', 2) can not be used here, which will change the polymerization type from free-radical polymerization to step-growth one.

(2) Step-growth polymerization

Formula: A + B -> A-B

Description: The function setMaxCris(string type0, int max_cris) will call step-growth polymerization model. The reaction could happen between active center A and reactive monomer B at the condition of both particles having h_cris =1 After reaction, the h_cris of B plus 1, the h_cris of A keeps it's value, and the h_cris of C minus 1.

1) Grafting to reaction

Formula: I* + ABBB -> I-ABBB

Description: The reaction happens between I and A.

reaction = galamost.Polymerization(all_info, neighbor_list, 2**(1.0/6.0), 16361)

reaction.setPr('I', 'A', 0.005)

reaction.setMaxCris('I', 1)

reaction.setMaxCris('A', 2)


app.add (reaction)

2) Polycondensation or polyaddition

Formula: A-B* + C-D -> A-B-C-D

reaction=galamost.Polymerization(all_info, neighbor_list, 2**(1.0/6.0), 123)

reaction.setPr('B', 'C', 0.005)

reaction.setMaxCris('B', 2)

reaction.setMaxCris('C', 2)


app.add (reaction)

(3) Exchange polymerization

Formula: A + B*-C->B*-A + C

Description: The function setExchangePr(string type0, string type1, string type2, float pr) will call exchange polymerization model. In a reaction, the pre-existed bond B-C will be disconnected and a new bond A-B will be generated simutanously.

1) Ligand exchange reaction

Formula: B + A*-C->A*-B + C

reaction = galamost.Polymerization(all_info, neighbor_list, 1.0, 123)

reaction. setExchangePr ('B', 'A', 'C', 0.1)

app.add (reaction)

Note: The A with h_init=1 is active center. B with h_cris=0 is reactive monomer. Before reaction, the h_cris of C shoud be h_cris>=1. After reaction, the h_cris of B plus 1, the h_cris of A keeps it’s value, and the h_cris of C minus 1.

2) Migration reaction

Formula: A + C*-A->C*-A + A

reaction = galamost.Polymerization (all_info, neighbor_list, 1.0, 123)

reaction. setExchangePr ('A', 'C', 'A', 0.1)

app.add (reaction)

Note: The particle C with h_init=1 is active center. Free particle A with h_cris=0 is reactive monomer. Before reaction, the h_cris of connected particle A shoud be h_cris>=1. After reaction, the h_cris of added A plus 1 and the h_cris of replaced A minus 1.


[1] Y.-H. Xue, Y.-L. Zhu, W. Quan, F.-H. Qu, C. Han, J.-T. Fan, H. Liu, Phys. Chem. Chem. Phys., 15, 15356, 2013.
[2] H. Liu, Y.-L. Zhu, J. Zhang, Z.-Y. Lu, Z.-Y. Sun, ACS Macro Lett., 1, 1249, 2012.
[3] H. Liu, Y.-L. Zhu, Z.-Y. Lu, F. Mueller-Plathe, J. Comput. Chem., 37, 2634, 2016.