Quantum molecular dynamic simulations of warm dense carbon monoxide

Quantum molecular dynamic simulations of warm dense carbon monoxide

Yujuan Zhang LCP, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China    Cong Wang LCP, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China    Dafang Li LCP, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China    Ping Zhang LCP, Institute of Applied Physics and Computational Mathematics, P.O. Box 8009, Beijing 100088, People’s Republic of China Center for Applied Physics and Technology, Peking University, Beijing 100871, People’s Republic of China

Using quantum molecular dynamic simulations, we have studied the thermophysical properties of warm dense carbon monoxide under extreme conditions. The principal Hugoniot, which is derived from the equation of state, shows excellent agreement with available experimental data up to 67 GPa. The chemical decomposition of carbon monoxide has been predicted at 8 GPa by means of pair correlation function. Based on Kubo-Greenwood formula, the dc electrical conductivity and the optical reflectivity are determined, and the nonmetal-metal transition for shock compressed carbon monoxide is observed around 43 GPa.

61.20.Ja, 51.70.+f, 31.15.A-, 64.30.-t
thanks: Corresponding author: zhang ping@iapcm.ac.cn

I Introduction

With rapid development of planetary science and explosion technology, pressure induced responses of materials under extreme conditions, where the combination of high temperatures and high pressures defines the warm dense regime, have recently attracted extensive attention Erns2009 (); Heml2000 (); Maze2004 (). Experimental techniques, e.g. gas gun Nell1981 (), chemical explosives Mint2006 (), pulsed power Knud2003 (); Boch2011 (), have been adopted in dynamic compressions, where pressures could reach megabar region. Due to the enormous progress in computational capacity, quantum molecular dynamics (QMD) Kress2004 (); Hols2008 (), provide powerful tools to study warm dense matter theoretically, where dissociations, recombinations, and ionizations characterize the high pressure behaviors.

As one of the major components in the reacted high explosives, the pressure introduced transitions of carbon monoxide (CO) are of crucial interest in exploring detonation products Kerl1980 (); Ross1980 (); Schm1997 (). The equation of states (EOS) and the relative quantities, for instance the Hugoniot points, sound velocity and particle velocity, are important parameters for materials under extreme conditions. Furthermore, electrical and optical properties such as, e.g., dc conductivity () and optical reflectance are closely related to the dielectric function, which can be evaluated through dynamic conductivity. The above mentioned quantities are of crucial interest in characterizing the unique behavior of warm dense CO. The Hugoniot points up to 70 GPa have been detected by a two-stage light-gas gun Nell1981 (), and three theoretical models were used to analyze the EOS and chemical decomposition. CO was reported to decomposes into condensed carbon and fluid CO at high pressures by pair potential calculations Ross1980 (), and the dissociation rates were then determined Alia2003 (). The principal Hugoniot and decomposition of CO have been investigated by means of QMD method Kress2004 (). However, wide range EOS, shock products, electrical and optical properties of CO under extreme conditions are still yet to be presented and understood.

In the present work, we have performed QMD simulations to study the thermophysical properties of warm dense CO. The decomposition and recombination of CO along the principal Hugoniot are analyzed by pair correlation functions (PCF). The dynamic conductivity is calculated by Kubo-Greenwood formula, from which the dc conductivity, the dielectric function and optical reflectivity can be extracted. The rest of the paper is arranged as follows. In section II, the computational method is briefly described. We present and discuss our calculated results, some of which are made comparison with available experimental and theoretical results in section III. Finally, we close our paper with a summary of our main results.

Ii Computational Method

In this study, the molecular dynamics trajectories are calculated by employing the Vienna simulation package (VASP) plane-wave pseudopotential code developed at the Technical University of Vienna Kres1993 (); Kres1996 (). The calculations are performed by a series of volume-fixed supercell with an invariable number of atoms. Employing the Born-Oppenheimer approximation, the electrons are fully quantum-mechanical treated by solving the Kohn-Sham equations for the orbitals and energies within a plane-wave, finite-temperature DFT formulation Leno2000 (), where the electronic states are populated according to the Fermi-Dirac distribution at electron temperature . The Perdew-Wang 91 (PW91) generalized gradient approximation Perd1991 () and the projector augmented wave (PAW) potential of Blöchl Bloc1994 () are employed to describe the exchange-correlation energy and the electron-ion interaction, respectively. The isokinetic ensemble (NVT) is employed for the ions, where the ionic temperature is controlled by Nośe thermostat Nose1984 () and the system is controlled in local thermodynamical equilibrium by setting the electron temperature and the ion temperature to be equal.

The supercell in our calculation contains 128 atoms (64 CO molecules) with periodic boundary condition. The plane-wave cutoff energy is selected to be 550 eV, which is tested to be good convergence for both total energy and pressure. In molecular dynamic simulation, only point of the Brillouin zone is included, while, 444 Monkhorst-Pack scheme Monk1976 () grid points are used in the electronic structure calculations. The selected densities range from 0.807 to 3.10 g/cm with temperatures from 77 to 11000 K, which highlight the principal shock Hugoniot region. All the dynamic simulations last 3-5 ps with time steps of 0.5-1.0 fs according to different conditions. The EOS data and pair correlation functions are obtained by averaging over all the particles and the final 1-2 ps simulations.

Iii results and discussion

iii.1 Equation of state and pair correlation function

(g/cm) (GPa) (K) (km/s) (km/s)
1.45 5.32 603 1.70 3.83
1.78 8.57 1614 2.40 4.39
2.25 20.20 4545 4.00 6.23
2.52 31.77 6600 5.17 7.60
2.70 43.27 8697 6.12 8.74
2.90 53.30 9175 6.90 9.56
3.10 67.27 11000 7.84 10.61
Table 1: Hugoniot pressure () and temperature () points derived from QMD simulations at a series of densities (). The corresponding particle velocity (), and shock velocity () are also displayed.

A precise description of EOS plays an important role in accurately calculating the electrical and optical properties. The EOS is examined theoretically along the Hugoniot, which can be derived from conservation of mass, momentum, and energy for a isolated system compressed by a pusher at a constant velocity. Rankine-Hugoniot equation describes the shock adiabat through a relation between the initial and final internal energies, pressures and volumes as follows Zeld1966 ():


where , and are the internal energy, volume, and pressure, and the subscripts 0 and 1 denote the initial and shocked states, respectively. In our canonical ensemble calculations, the internal energy equals to the sum of the total energy from the finite-temperature DFT calculation and zero-point energy. The pressure consists of the electronic and ionic components. The pressure can be obtained by , where is the ion number density, and is the Boltzman constant.

Figure 1: (Color online) Principal Hugoniot (the main panel) and the (, ) diagram (the inset) of CO. For comparison, the present QMD results, together with the previous experimental data (Ref. 4) and other theoretical results (Ref. 4 and Ref. 12) are all provided.

As the starting point along the Hugoniot, the initial density is =0.807 g/cm at a temperature of 77 K Kress2004 (), with the relative internal energy eV/molecule. The initial pressure can be neglected compared to the high pressure of the final state. The Hugoniot points are determined in the following way. For a given density , a series of simulations are executed for different temperature . and as a function of temperature are then fitted to polynomial expansions of . The principal Hugoniot temperature and pressure can be determined by solving Eq. (1). The particle velocity and shock velocity are then determined from the other two Rankine-Hugoniot equations Zeld1966 (), and . The principal Hugoniot points of CO derived from Eq. (1) as well as the particle velocity and shock velocity are listed in Table I.

The simulated principal Hugoniot points for CO are shown in Fig. 1, where previous experimental and theoretical data are also provided for comparison. As seen in Fig. 1, our simulated results accord excellently with the experimental results Nell1981 () and the theoretical results Kress2004 (). The (, ) diagram also shows good agreement with experimental results (the inset of Fig. 1).

Figure 2: (Color online) Pair-correlation functions for C–C (black line), O–O (red line), C–O (blue line) along the principal Hugoniot of CO.

To examine the structural transitions of CO under shock compressions, the PCFs [=], which represent the possibility of finding a particle at a distance from a reference atom, are calculated. The PCFs and the corresponding atomic structures along the principal Hugoniot of CO, are presented in Fig. 2 and Fig. 3, respectively. At the pressure of 5 GPa (=1.45 g/cm and =603 K), the peak of C–O PCF g() locates around 1.13 Å, which corresponds to the equilibrium internuclear distance of the C–O bond in CO molecule Gill1950 () [Fig. 2(a)], indicating that CO molecules remain their ideal molecular fluid state without dissociation. From the corresponding atomic structure [Fig. 3(a)], it can be clearly seen that the supercell is filled with CO molecules. With the increase of pressure along the Hugoniot, the main peak of PCF g() becomes reduced in amplitude and gets broadened due to the thermal excitation, as shown in Fig. 2(b). Meanwhile, the g() appears a peak at =1.47Å, which lies between typical and hybrid C–C bond length. From the atomic structure of CO at =2.25 g/cm (=20.2 Gpa, =4545 K), CO molecules and the small clusters with carbon backbones can be clearly observed [Fig. 3(b)]. The structural transitions of CO may bring with the increase of the , which will be discussed in the following section. As the difference between the length of C-O bond in CO molecule (1.13Å) and that of CO molecule (1.16Å) wang2010 () is very small, the appearances of their PCFs are similar. Despite the fact that the formation of CO increases the peak value of g(), CO is further dissociated leading to a decrease in the peak value. Under the competition of the above two effects, the amplitude of PCF g() still keeps reduced as the pressure increases. At higher pressures, CO molecules get dissociated, and a mixture of CO, atomic carbon, and some carbon backbones clusters forms, which results in a further decrease of the peak value of PCF g(). In particular, diatomic oxygen [the peak of PCF g() is around 1.23 Å as shown in Fig. 2(d)] and monoatomic oxygen are formed in the shocked system [Fig. 3(d)], which lead to the nonmetal-metal transition. Overall, the systematic behavior of the dissociation and recombination of fluid CO under shock pressure is described. At low pressures (20 GPa), CO as well as some carbon backbone clusters forms. With increasing the pressure (43GPa), CO dissociates and a mixture of CO, atomic carbon diatomic oxygen monoatomic oxygen and carbon backbone clusters forms.

Figure 3: (Color online) Atomic structures at different Hugoniot points. The carbon and oxygen atoms are denoted by gray and red balls, respectively, and the relative isosurface of charge density (blue regimes) are plotted. (a)=1.45 g/cm, =603 K; (b)=2.25 g/cm, =4545 K; (c)=2.52 g/cm, =6600 K; (d)=2.70 g/cm, =8697 K.

iii.2 Dynamic conductivity under shock compressions

The dissociation of molecules under dynamic compression are closely related to the nonmetal-metal transition properties. The dynamic conductivity is derived from the Kubo-Greenwood formula as follows Kubo1957 (); Gree1958 ():


where is Kohn-Sham eigenstate, with corresponding eigenvalue , and occupation number . is the volume of the supercell. (k) is the K-point weighting factor. The and summations range over discrete bands included in the calculation, and over the three spatial directions.

In order to calculate the conductivity, thirty atomic configurations are selected at equilibrium along the principal Hugoniot. The behavior of the conductivity versus energy at different Hugoniot points are calculated (not shown). With the increase of density and temperature along the principal Hugoniot, the main peak increases in amplitude and gets broadened. Such change of the peak indicates the increase of another important physical quantity, dc conductivity, which follows the static limit . The dc conductivity is extracted as shown in Fig. 4. For the pressures below 8 GPa, can be treated as small as negligible, which shows the insulating nature of CO at low pressures. At higher pressures, the dc conductivity of CO jumps three orders of magnitude from 8 to 67 GPa (603 to 11000 K). The nonmetal-metal transition occurs around 43 Gpa, which is accompanied by the gradual dissociation of the molecular fluid. Such nonmetal-metal transition under shock pressure has also been reported in fluid CO wang2010 () and O Wangc2010 ().

Figure 4: The dc conductivity dependence on pressure is extracted from the dynamic electrical conductivity.

The imaginary part of conductivity can be obtained via the Kramers-Kronig relation as


where is the principal value of the integral and is frequency. The real part () and the imaginary part () of dielectric functions can be derived from the two parts of the conductivity,




The square of the index of refraction, which contains the real part () and the imaginary part (), is equal to the dielectric function. Then the index of refraction can be obtained by



Figure 5: (Color online) Optical reflectivity of shocked CO for wavelengths of 404, 808 and 1064 nm along the principal Hugoniot.

The index of refraction is useful for evaluating optical properties, such as the reflectivity and absorption coefficient. Emissivity of spectrum is an important quantity, which can be used to determine temperature in experiments. The optical reflectivity and emissivity of spectrum are closely related. We calculated the reflectivity of CO under shock pressure by the following relation


The optical reflectivities of shocked CO for wavelengths of 404, 808 and 1064 nm along the principal Hugoniot are shown in Fig. 5. The change induced by the pressure can be clearly seen. The optical reflectance increases from 0.12 to 0.46-0.60 with the pressure increasing from 8 to 67 GPa, which can be interpreted as a gradual transition from nonmetal state to metal state. Similar behavior in the optical reflectance of liquid deuterium has been previously observed in experiment Cell2000 (). Our predictions of CO need be tested in the future experiments.

Iv Conclusions

In summary, we have performed QMD simulations to study the thermophysical and optical properties of shock-compressed carbon monoxide along the principal Hugoniot. The principal Hugoniot curve has been derived from the EOS, which shows excellent agreement with experimental and theoretical results. The dissociation of CO under shock compressions is described by the pair correlation functions. At low pressures (20 GPa), CO and large carbon backbone clusters form. Above 43 GPa, CO molecules dissociate and a mixture of CO, atomic carbon and carbon backbones clusters form. Based on the calculation of the dynamic conductivity, the nonmetal-metal transition and the change in optical reflectivity are observed. In particular, the nonmetal-metal transition occurs around 43 GPa.

This work was supported by NSFC under Grants No. 11005012 and No. 51071032, by the National Basic Security Research Program of China, and by the National High-Tech ICF Committee of China.


  • (1) R. Ernstorfer, M. Harb, C.T. Hebeisen, T. Dartigalongue, and R.J.D. Miller, Science 323, 1033 (2009).
  • (2) R.J. Hemley, Annu. Phys. Chem. 51, 763 (2000).
  • (3) S. Mazevet, P. Blottiau, J.D. Kress, and L.A. Collins, Phys. Rev. B 69, 224207 (2004).
  • (4) W.J. Nellis, F.H. Ree, M. van Thiel, and A.C. Mitchell, J. Chem. Phys. 75, 3055 (1981).
  • (5) V.B. Mintsev and V.E. Fortov, J. Phys. A 39, 4319 (2006).
  • (6) M.D. Knudson, D.L. Hanson, J.E. Bailey, C.A. Hall, J.R. Asay, Phys. Rev. Lett. 90, 035505 (2003).
  • (7) I.A. Bocharova, A.S. Alnaser, U. Thumm, T. Niederhausen, D. Ray, C.L. Cocke, and I.V. Litvinyuk, Phys. Rev. A 83, 013417 (2011).
  • (8) J.D. Kress, S. Mazevet, L.A. Collins, and P. Blottiau, AIP Conference Proceedings, 706, 289 (2004).
  • (9) B. Holst, R. Redmer, and M. P. Desjarlais, Phys. Rev. B 77, 184201 (2008).
  • (10) G.I. Kerley and J.J. Abdallah, J. Chem. Phys. 73 5337 (1980).
  • (11) M. Ross and F.H. Ree, J. Chem. Phys. 73, 6146 (1980).
  • (12) S.C. Schmidt, D.S. Moore and M.S. Shaw, J. Chem. Phys. 102 325 (1997).
  • (13) A. Aliat, A. Chikhaoui, and E.V. Kustova, Phys. Rev. E 68, 056306 (2003).
  • (14) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  • (15) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  • (16) T. Lenosky, S. Bickham, J. Kress, and L. Collins, Phys. Rev. B 61, 1 (2000).
  • (17) J.P. Perdew, Electronic Structure of Solids (Akademie, Berlin, 1991).
  • (18) P.E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  • (19) S. Nosé, J. Chem. Phys. 81, 511 (1984).
  • (20) H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13 5188 (1976).
  • (21) Y. Zeldovich and Y. Raizer, Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena (Academic Press, New York, 1966).
  • (22) O.R. Gilliam, C.M. Johnson, and W. Gordy, Phys. Rev. 78, 140 (1950).
  • (23) C. Wang and P. Zhang, J. Chem. Phys. 133, 134503 (2010).
  • (24) R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957).
  • (25) D.A. Greenwood, Proc. Phys. Soc. London 71, 585 (1958).
  • (26) C. Wang and P. Zhang, J. Chem. Phys. 132, 154307 (2010).
  • (27) P.M. Celliers, G.W. Collins, L.B. Da Silva, D.M. Gold, R. Cauble, R.J. Wallace, M.E. Foord, and B.A. Hammel, Phys. Rev. Lett. 84, 5564 (2000).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description