A Study of thermal decomposition in cellulose by molecular dynamics simulation

PM3 method was used in this paper to optimize cellulose molecular structure which is the main component of biomass and a series of structural parameter was attained. The single chain of cellulose (the degree of polymerization is 9) was simulated in different force fields by molecular dynamic method. Energy history, deposition temperature and the cracked groups of simulation process in different force fields was gotten, of which Amber force field is quite matched to the experiments data. By simulating the process of cellulose thermal decomposition with MD which is based on Amber force field and quantum mechanics, we get the sequence of bond break of cellulose molecule and the first cracked group. Also, the first production was analyzed. The heating process includes two stages: vibrate at low temperature and break at high temperature (273k-375k) and breaking stage when the temperature of the system arrived at 375K.


INTRODUCTION
The special meeting, green energy: the cooperation of government and enterprise, of Boao Asian Forum pointed out that the miraculous development of Asia will not continue if no appropriate measures were taken to make sure energy safety, decrease energy consuming, find new energy and to release environment stress.However, all the problems probably can be solved by biomass energy, which has become one of the most popular and important topics in the word.
Biomass consists of cellulose, hemicellulose, lignin and a small amount of ash content.Cellulose, which is D-glucose high molecular polymer formed by connection between  (1-4)-glycosidic bond, is the main component of biomass, accounting for 40%-96% of the total amount of biomass.Many researches on relation between raw material and production of biomass energy were done by scholars form all over the world, using thermal decomposition, liquefaction, gasification [2,3,4,5].However, the work on the process is rare.In order to investigate the further principle of biomass energy's thermal deposition process, a reactive molecular dynamics model was developed to simulation the thermal deposition of cellulose in this paper.
Molecular structures were optimized before molecular dynamic simulation to lower the molecule energy to the possibly lowest degree.The lower the energy is, the steadier the structure is and the higher probability exists in the system is.Single point calculation was applied for optimizing energy; Optimized structure was searched by calculating a series of bond length and bond angle.
Single molecular was employed in this MD simulation, which differs form other poly-molecular systems.The lowest point of partial total potential energy was gotten when searching the lowest total potential energy, but the global optimization was needed.MD method must be used when optimal three dimensional structures were set.Molecules are heated up and the structure extends and relaxes adequately at high temperature.Then we cool down the molecules and calculate the optimum structure.Optimal three dimensional structures can be obtained.The conformations were much different between single molecule and multiple molecules which are more accordant with practical situation.Single molecule was researched in this paper considering huge system, complex structure and the operation ability of computes.
The semi-empirical method MNDO-PM3 based on MNDO model was used in this paper and Polak-Ribiere conjugate gradient was applied in optimization with RMS setted as 0.042kJ/mol.The results of cellulose molecular optimization were listed in Table 1.

Simulation Model
The molecular formula of cellulose is (C6H10O5) n , where n is polymerization degree.Haworth structure is used to express the cellulose molecule in the simulation [6].The structure model is given in Fig. 1.
The simulated object is single chain of cellulose with polymerization degree is 9.The original size is 6.9727 ).Periodic boundary condition was applied.During the simulation, the temperature was heated up from the 293K at the begging to 1273K which is simulation temperature.The heating time is 100ps; the simulation time is 10ps; the step size is 0.001ps.The parameters of energy and temperature were collected every time step.

Assumptions in Simulation
1) The broken groups have no influence on subsequent bond break.
2) The broken groups didn't get together and form new molecules.

FOFCE FIELD
Force field, with relation to reliability of the result, is the base of molecular dynamic simulation.The force fields become more and more complex as computing systems swell.It has different forms with its advantages and limits in each form.Some force fields can be chosen when the organic molecule was simulated, such as MM + , AMBER, CHARMM, OPLS [7,8,9,10].

Force Fields for Simulation
The total potential energy expression of the macromolecule is partitioned into several energy terms as given in Eq.1.These contributions include non-bonding energies (U nb ), bond stretching energies (U b ), angle bending energies, ( U  ), torsion angle energies U  , out-of-plant bending energies ( x U ), columbic interaction energies ( el U ).
AMBER (Assisted Model Building and Energy Refinement) force field which was developed by Peter A Kollman and coworkers in University of California San Francisco was widely used for proteins and DNA.Force field functions and parameter sets are derived from both experimental work and high-level quantum mechanical calculations.The first three terms in Eq.2 denote the internal coordinates of bond stretching, angle bending and torsions.The non-bonded terms account for the Van DerWaals and electrostatic interactions and the last term represent the 12-6 Lennard-Jones hydrogen bond treatment.
(2) CHARMM (Chemistry at Harvard Macromolecular Mechanics) force field was developed by Harvard, the parameter of which come form not only the result of experiment and computing, but also calculation result of quantum.It was used for many molecular systems, including organic micro molecule, solution, polymer, biochemical molecule.The form of the potential energy function we will use is given by the following equation [3]: ( , , ) ij on off sw r r r is switching function.
MM + force field was developed by Allinger.El, some common atoms are divided different forms, which have different parameters.It was applied in organic compounds, free radical, ions.The normal form is given: is cross action term which is given by following EQ while the bond lengths are 1 r and 2 r for the same molecule.OPLS force field formed by OPLS-UA model and OPLS-AA model was used for simulating organic molecules and multi-peptide.The parameters of bond extension and bend were attained form AMBER force field.OPLS force field is mainly used for computing the conformation energy of gaseous organic molecule, hydration free energy of organic liquid, and other thermodynamics features.

The Energy History
During the heating process (100ps＞t＞0ps), the total energy increases as temperature rises.But in the simulation process (t＞100ps), the temperature is constant and the total energy is steady.Form Fig. 2, as can be seen, the total energy lowered to minimum in a short time at the beginning of simulation.The energy was adjusted after the system was optimized and before simulating, the range of which was E o ＞E c ＞E m ＞E a .Table 2.
shows the energy adjustment range of OPLS is largest and the total energy is E ch ＞E am ＞E mm ＞E op .

Comparison of Different Force Fields
The straight chain starts to bend from both sides to the middle as the energy of each atom in the chain of cellulose rises; atoms move and vibrate more and more strongly.The groups in the cellulose begin to break when the total energy of system come to the special value.The breaking temperature is TO＞TC ＞TA＞TM as show in Table 3. Cellulose breaks more and more strongly as the temperature rises while only a few groups break at the beginning.The temperature of main decomposition process is shown in Table 2, from which the conclusion can be drawn: AMBER force field quite agrees with the experiment value [11].

The Energy History in Heating Process
The system was heated at the initial temperature of 273K.Fig. 3 shows the whole temperature history during the heating-up process.The temperature gradually rises during heating process with bigger fluctuation at higher temperature.The total energy history is shown in shown in Fig. 4, form which the rise of total energy as time going can be seen.

Breaking of Molecular Chain
Not considering bond bonding, the heating process includes two stages: vibrate at low temperature and break at high temperature.During low temperature (273k-375k), the length of bond grows, bond angles bend, torsion angles increase.Stretching energies (U b ), angle bending energies (U θ ), torsion angle energies(U Φ ) and out-of-plant bending energies (U χ ) increase, but the bond was not broken due to non-bonding energies (U nb ) and columbic interaction energies (U el ).As seen in Fig. 5, atoms vibrate more and more strongly and   the straight chain starts to bend from both sides to the middle.When the temperature of the system arrived at 375K, it comes to breaking stage.The groups move stronger with the increasing of temperature, some of which overcome the electronic force and van der waals force and leave from the long chain.Then, chemical bonds break and cellulose starts to decompose (Fig. 6).
When the length of chemical bond is larger than 1.2 times of its original length, we think that it is broken or disappears.Many groups are formed when cellulose begins to break and many of short chain groups can be obtained in low temperature process, such as: A lot of free groups [O] accelerate the oxidation reaction.Bottom (-OH) is oxidated to corresponding aldehyde and acids.For example: The re-combination of these broken groups form the first production such as: CO 2 , H 2 O(L), CH 4 , alkanes such as:CH 3 -CH 3 , olefin such as CH 2 =CH 2 , aldehydes such as CH 2 O.With the further increasing of temperature, organic biological oil with more than 6 carbons formed with highest production rate at around 800K-850K, which is quite matched with the experiments [11].

The Order of Molecule Breaking
The order of cellulose bonds breaking are obtained when the broken groups are shielded after the cellulose breaking.Fig. 7 shows the planar structure of carbon and oxygen, and the hydrogen is not shown because the breaking of (O-H) and (C-H) is not involved.The numbers show the orders and the letters show the units number of cellulose.
In one unite of cellulose, the hydroxyl groups (-OH) in the ring shed first, then the hydroxyl (-OH) in branched chain and the ring.Sometimes the whole ring breaks directly, unit I, for example.The (C-O) which has the lowest energy in the chain of cellulose breaks first because decomposition always occurs from the lowest energy.
As for the whole chain of cellulose molecule, the molecule chain decomposes from both sides to the middle gradually.The hydroxyl (-OH) of inside unit will break earlier than the ring of two-terminals.Separate adjustment is reasonable because the order given in this paper is not steady for the randomness of chemical reaction.In spite of thus bugs, the general tendency is correct.

CONCLUSIONS
1) Series of parameters of cellulose structure were obtained by optimizing the cellulose chain.
2) AMBER force field is more adaptive for simulating cellulose with lower polymerization degree than others.
3) Details of cellulose thermal decomposition and the ranges of decomposition temperature are gotten from the simulation.The order of cellulose unit breaking is obtained.Therefore, the detailed process of cellulose thermal decomposition is shown.
4) The heating process includes two stages: vibrate at low temperature and break at high temperature (273k-375k) and breaking stage when the temperature of the system arrived at 375K.
5) The re-combination of these broken groups which were gotten from thermal decomposition forms the first production which was the most important product.
The conclusions concluded from the simulation of cellulose thermal decomposition by molecular dynamic model matches the experiment quite well, convincing molecular dynamics could be a very significant tool for science research.The results of the simulation are very significative for the following researches.

Figure 2 .
Figure 2. The energy histories of simulation process.

Figure 3 .
Figure 3.The temperature history of system.Figure 4. The energy history of system.

Figure 4 .
Figure 3.The temperature history of system.Figure 4. The energy history of system.

2
) in the middle also can be oxidated to ketone or be alkene off-(-OH),

Figure 5 .
Figure 5.The process of heating.Figure 6.The process of decomposition.

Figure 6 .
Figure 5.The process of heating.Figure 6.The process of decomposition.

Figure 7 .
Figure 7.The breaking order of cellulose single chain.

Table 1 .
The parameters of cellulose molecule before and after optimizing.

Table 2 .
Different Force Fields Parameters.