A Study of thermal decomposition in cellulose by molecular dynamics simulation

doi:10.4236/ns.2009.11008

Paper Menu >>

Journal Menu >>

Natural Science, 2009, 1, 41-46 NS

http://dx.doi.org/10.4236/ns.2009.11008

A Study of thermal decomposition in cellulose by

molecular dynamics simulation

Chao Liu1, De-Zheng Jiang1, Shun-An Wei2, Jin-Bao Huang1

1College of Power Engineering, Chongqing University, Chongqing 400044, China; 2College of Chemistry and Chemical Engineering,

Chongqing University, Chongqing 400044, China.

Email: [email protected]

Received 19 May 2009; revised 5 June 2009; accepted 7 June 2009.

ABSTRACT

PM3 method was used in this paper to optimize

cellulose molecular structure which is the main

component of biomass and a series of struc-

tural parameter was attained. The single chain

of cellulose (the degree of polymerization is 9)

was simulated in different force fields by mo-

lecular dynamic method. Energy history, depo-

sition 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 ar-

rived at 375K.

Keywords: Molecular Dynamics Simulation; Cel-

lulose; Thermal Decomposition; Bond Breaks

1. 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 connec-

tion 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 dynam-

ics 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 simula-

tion, 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 struc-

tures were set. Molecules are heated up and the struc-

ture extends and relaxes adequately at high tempera-

ture. 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.

42 C. Liu et al. / Natural Science 1 (2009) 41-46

2. SIMULATION TECHNIQUE

2.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

(





)×6.2610(





)×41.3464( 



) after optimizing. The

size of the simulation cell is 15(



)×15(



)×50(





Periodic boundary condition was applied. During the

simulation, the temperature was heated up from the

293K at the begging to 1273K which is simulation tem-

perature. The heating time is 100ps; the simulation time

is 10ps; the step size is 0.001ps. The parameters of en-

ergy and temperature were collected every time step.

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

Parameters

Groups Before Optimizing After Optimizing Parameters

Groups Before Optimizing After Optimizing

C1H1 1.09A 1.1073A O2H 0.96A 0.9617A

C1H2 1.09A 1.10458A C3C4C5 108.247° 110.169°

H1C1H2 109.47° 108.014° C4O2H 109.471° 108.112°

C1O1 1.43A 1.40063A HC4O2H -159.4365° -159.246°

O1H1 0.96 A 0.96587A C3C4C5O2 -119.811° -127.222°

C1C2 1.54A 1.55006A C3C4HC5 119.8° 120.336°

HO1C1 109.471° 107.203° C5H 1.09A 1.11574A

HC1HO1 -120° -116.391° C5O3 1.43A 1.40229A

C2H 1.09A 1.11871A O3H 0.96A 0.95074A

C2O4 1.4326A 1.42258A C5C6 1.54012A 1.55571A

C2C3 1.54A 1.54324A C4C5C6 108.867° 112.894°

C3C2H 109.325° 109.75° C5O3H 109.471° 107.744°

C3C2O4 110.057° 115.958° C4C6C5H -120.01° -119.849°

HC3C2O4 120.044° 111.173° C4C6C5O3 119.99° 128.836°

HC3C2C1 119.857° 121.745° C6H 1.09A 1.12113A

C3H 1.09A 1.1336A C6O5 1.43A 1.41851A

C3C4 1.53748A 1.5499A C6O4 1.43285A 1.39737A

HC3C2 109.62° 109.505° C5C6O4 110.052° 116.06°

C2C3C4 108.875° 113.707° C5C6O5 109.339° 113.56°

HC3C4 109.62° 109.879° C5C6O5H 119.728° 127.655°

C4HC2C3 38.1592° 34.3572 C5C6O4O5 120.104° 123.063°

C4O2 1.43A 1.40198A C5C6HO5 119.743° 127.236°

C4H 1.09A 1.12203A C6O4C2 119.743° 120.239°

C4C5 1.53753A 1.54967A

Figure 1. The model of cellulose molecule (n=9).

C. Liu et al. / Natural Science 1 (2009) 41-46 43

2.2. 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.

3) The broken groups didn’t decompose secondarily.

3. 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 lim-

its 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].

3.1. Force Fields for Simulation

The total potential energy expression of the macromole-

cule is partitioned into several energy terms as given in

Eq.1. These contributions include non-bonding energies

(Unb), bond stretching energies (Ub), angle bending en-

ergies, (U



), torsion angle energies U



, out-of-plant

bending energies (

U), columbic interaction energies

(el

U).

nb bel

UUU UUUU



 (1)

AMBER (Assisted Model Building and Energy Re-

finement) 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 treat-

ment.

0000

()( )[1cos()]

UKbb KVn











 

12612 10

[( )2( )]()

ij ijij ij

ij ij

rrrrCrD r









(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, in-

cluding organic micro molecule, solution, polymer, bio-

chemical molecule. The form of the potential energy

function we will use is given by the following equation

[3]:

()( )[||cos()]

Ukrr kkkn



 









22 2

012 6

()( )(,,)

i jijij

ijon off

ij ij

ij ij ij

qqA B

kswrrr

rrr



 

 



(3)

where 22 2

(, ,)

ijon off

wr rr 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 com-

pounds, free radical, ions. The normal form is given:

nb belcross

UUU UU UUU





 

(4)

where 6

()( )

Ur aebr







 

()(1 cos)







 ()(1 cos2)

Ux kx

223

0000

()()[1()()()]

Urkkk





 

 

 

223

0000

()()[1()() ()]

bbbb

Urrrkrrkrrkrr



 

cross

U is cross action term which is given by follow-

ing EQ while the bond lengths are 1

r and 2

r for the

same molecule.

1 211,022,0

(, )()()

Urrr rrr

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, hy-

dration free energy of organic liquid, and other thermo-

dynamics features.

()()

[1 cos()][1 cos(2)]

reqeq

Ukrr kV









 

 



2126

[1 cos(3)]()

i jijijijijij

qq erArCr



 



(5)

3.2. The Energy History

During the heating process (100ps＞t＞0ps), the total

energy increases as temperature rises. But in the simula-

tion 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

44 C. Liu et al. / Natural Science 1 (2009) 41-46

after the system was optimized and before simulating,

the range of which was Eo＞Ec＞Em＞Ea. Table 2.

shows the energy adjustment range of OPLS is largest

and the total energy is Ech＞Eam＞Emm＞Eop.

3.3. 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].

4. SIMULATION RESULTS BASE ON

AMBER FORCE FIELD

4.1. 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

Figure 2. The energy histories of simulation process.

in Fig. 4, form which the rise of total energy as time

going can be seen.

4.2. 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 (Ub),

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

(Unb) and columbic interaction energies (Uel). As seen

in Fig. 5, atoms vibrate more and more strongly and

Table 2. Different Force Fields Parameters.

AmberCharmmMM+ Opls Experiment

Initial Energy (kJ/mol) 2644.3 3606.6 1924.6 4418.3

Adjusted Energy (kJ/mol) 2502 3117.1 2288.6 1912.1

Energy difference (kJ/mol) -142.3 -489.5 364 -2506.2

Balanced Energy (kJ/mol) 7564.7 7836.6 7455.9 6757.2

Temperature of starting decomposition (K)380 390 310 417 400

Main Temperature Ranges of

Decomposition (K) 420~750450~860310~710600~950 480~700

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

C. Liu et al. / Natural Science 1 (2009) 41-46 45

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: (-OH), (-CO-), (-CHOH-

CHOH-), (-CHOH-CH2-CH2-), (-CH2OH), (-CH3),

(-CH2-CH2-), (-CH2-CH2-CHOH-), (-CH=CH-CHOH-),

(-CH=CH-), (-CHOH-CHOH-CH3), (-CH2-CH2-CH2- CH2-),

(-CHOH-CH2-CH(C)-CHOH-), (-CHOH-CH2-CH2-). A

number of unsteady (OH) groups which are prone to re-

lease free groups [O] exist in the system after breaking.

()()[]OHOHH OO  (6)

A lot of free groups [O] accelerate the oxidation reac-

tion. Bottom (-OH) is oxidated to corresponding alde-

hyde and acids. For example:

[]

323 2

CH CHOHCH CHOHO

(7)

[]

CH CHOCH COOH 

(8)

(-OH) in the middle also can be oxidated to ketone or be

alkene off-(-OH),

[]

33332

() ()

CHCOHHCHCHCOCHH O

(9)

322 22

CHCHOHCHCHH O (10)

and:

[]

CO CO (11)

33 33

CH CHCH CH



 (12)

CHOHCHOH



 (13)

[]

422

CHCH OH (14)

The re-combination of these broken groups form the

first production such as: CO2, H2O(L), CH4, alkanes

such as:CH3-CH3, olefin such as CH2=CH2, aldehydes

such as CH2O. With the further increasing of tempera-

ture, organic biological oil with more than 6 carbons

formed with highest production rate at around

800K-850K, which is quite matched with the experi-

ments [11].

4.3. The Order of Molecule Breaking

The order of cellulose bonds breaking are obtained when

the broken groups are shielded after the cellulose break-

ing. Fig. 7 shows the planar structure of carbon and oxy-

gen, 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 cel-

lulose.

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.

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

46 C. Liu et al. / Natural Science 1 (2009) 41-46

Figure 7. The breaking order of cellulose single chain.

5. CONCLUSIONS

1) Series of parameters of cellulose structure were ob-

tained 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 ob-

tained. Therefore, the detailed process of cellulose ther-

mal 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 cel-

lulose 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.

ACKNOWLEDGEMENTS

The Project is sponsored by the National Natural Science Foundation

of China (No. 50776101)

REFERENCES

[1] Cook, J. and Beyea, J. (2000) Bioenergy in the United

States: Progress and possibilities. Biomass and Bio-

energy, 18(6), 441-441.

[2] McKendry, P. (2002) Energy production from biomass:

gasification technologies. Bioresource Technology, 83,

55-63.

[3] McKendry, P. (2002) Energy production from biomass:

Conversion technologies. 83, 47-54.

[4] Bridgwater, A. V. (1999) Principles and practice of bio-

mass fast pyrolysis processes for liquids. Journal of Ana-

lytical and Applied Pyorlysis, 51, 3-22.

[5] Slesser, M. and Chris, L. (1982) Biological energy re-

sources. London: E&F. N. Spom Ltd.

[6] Purves, C. B. (1954) Cellulose and cellulose derivatives.

Part 1, Interscience, New York, 29-98.

[7] Weiner, S. J. and Kollman, P. A. (1984) D A Case. J. Am.

Chem. Soc, 106, 765-784.

[8] Brooks, B. R., Bruccoleri, R. E., Olafson, B. D. et al.

(1983) J. Comput. Chem, 4, 187-217.

[9] Allinger, N. L. and Yan, L. (1993) J. Am. Chem. Soc.,

115-119.

[10] Qui~nonero, D., Tomas, S. and Frontera, A. (2001)

OPLS all-atom force field for squaramides and squarica-

cid. Chemical Physics Letters, 9, 331-338.

[11] Liao, Y. F. (2003) Mechanism study of biomass pyroly-

sis. [Ph. D. Thesis]. Hang Zhou, Zhejiang University,

China.

[12] Tan, H. (2005) Mechanism study of biomass pyrolysis.

[Ph. D. Thesis]. Hang Zhou, Zhejiang University, China.