氢氘物态方程研究进展

刘海风 张弓木 张其黎 宋红州 李琼 赵艳红 孙博 宋海峰

引用本文:
Citation:

氢氘物态方程研究进展

    作者简介: 刘海风(1968—), 女, 硕士, 研究员, 长期从事物态方程理论研究. E-mail: liu_haifeng@iapcm.ac.cn;
  • 基金项目: 国家自然科学基金青年项目 11604018
    中国工程物理研究院科学技术基金 990105
    国家自然科学基金重点项目 10135010
    科学挑战计划 TZ2016001

  • 中图分类号: O521.2

Progress on Equation of State of Hydrogen and Deuterium

  • CLC number: O521.2

  • 摘要: 针对近二十多年的氢氘物态方程理论研究工作进行了综述分析,结合本课题组改进的自由能模型、直接量子蒙特卡洛和量子分子动力学方法的模拟结果,对多个研究小组采用不同方法获得的氢氘宽区物态方程数据进行了定量评估分析。结果表明:在当前理论框架下,仅基于第一原理数值模拟得到的氢氘物态方程能够描述的热力学相空间有限;多模型集成的H-REOS.3数据库在105 K以下温度与数值模拟结果的相对差别较大,且数据稀少,二者均不能满足工程应用需求。建议采用基于半经验模型的宽区域物态方程研究方法,即结合高精度的实验研究、数值模拟和解析模型,构建满足工程应用需求的氢氘宽区实用物态方程。
  • 图 1  QMD模拟氘物态方程的压强相对偏差

    Figure 1.  Relative differences of pressure from QMD simulation for deuterium

    图 2  QMD模拟氢物态方程的压强相对偏差

    Figure 2.  Relative differences of pressure from QMD simulation for hydrogen

    图 3  不同泛函计算的氘物态方程压强相对差别

    Figure 3.  Relative differences of pressure from QMD simulation for deuterium in different exchange-correlation analysis

    图 4  PIMC模拟氘物态方程的压强相对偏差

    Figure 4.  Relative differences of pressure from PIMC simulation for deuterium

    图 5  PIMC模拟氢物态方程的压强相对偏差

    Figure 5.  Relative differences of pressure from PIMC simulation for hydrogen

    图 6  氢的MP2011与REOS物态方程的相对偏差

    Figure 6.  Relative differences of equation of state for hydrogen between MP2011 and REOS

    图 7  RPIMC模拟的氘物态方程与REOS物态方程的相对偏差

    Figure 7.  Relative differences of equation of state for deuterium between RPIMC and REOS

    图 8  DFT模拟的氢物态方程与REOS物态方程的相对偏差

    Figure 8.  Relative differences of equation of state for hydrogen between DFT and REOS

    图 9  改进的自由能模型计算的氢物态方程(LiEOS)与REOS物态方程的相对偏差

    Figure 9.  Relative differences of equation of state for deuterium between modified free energy model (LiEOS) and REOS

    表 1  氢氘物态方程的QMD模拟工作

    Table 1.  Works of QMD simulation for equation of state of hydrogen and deuterium

    Isotope Time Firstauthor Ref. Method Code Organization Temperature/K Density/(g·cm-3)
    D 1997 T. J. Lenosky [135] TBMD LANL 3 000-31 250 0.58-1.47
    D 2000 G. Galli [136] CPMD GP code LLNL 10 000 0.67,1.0
    D 2000 S. Bagnier [137] BOMD VASP CEA The effect of the local-spin-density-approximation functional
    D 2000 T. J. Lenosky [138] BOMD VASP LLNL,LANL 2 000-31 500 0.506-0.851
    D 2002 F. Gygi [139] CPMD GP code LLNL The compressibility is determined by shock-induced electronic excitations
    D 2003 M. P. Desjarlais [140] BOMD VASP SNL 3 800-50 105 0.553-1.756
    D 2004 S. A. Bonev [141] CPMD GP code LLNL 20-19 860 0.171-0.779
    D 2004 H. F.Liu [149] BOMD VASP IAPCM 1 000-10 000 0.506-1.0
    D 2016 J. F. Danel [142] QMD,OFWMD ABINIT,VAAQP CEA 11 604-116 045 0.2-20
    D 2016 V. V. Karasiev [143] QMD,OFMD PROFESS@Q-ESPRESSO,ABINIT Universityof Florida 2 000-250 000 0.2-10, No data list
    D 2017 M. D. Knudson [12] QMD VASP,Diff Exc SNL, WashingtonState University
    H 2001 L. A. Collins [146] BOMD VASP LANL,LLNL 5 000-30 000 0.334-0.525
    H 2008 B. Holst [147] BOMD VASP Institut für Physik,SNL 500-20 000 0.5-5.0
    H 2010 M. A. Morales [126] BOMD,QMC QBOX,CEIMC University of Illinois,University of L'Aquila 2 000-10 000 0.724-2.329
    H 2011 L. Caillabet [148] QMD,CEIMC PIMC, QMD, CEIMC CEA -116 045 0.2-5
    H 2013 C. Wang [151] QMD,OFMD ABINIT IAPCM 1.564×104-5.004441×107 9.82×10-3-1.347×103
    下载: 导出CSV

    表 2  氢氘物态方程的QMC数值模拟

    Table 2.  Works of QMC simulation for equation of state of hydrogen and deuterium

    Isotope Time Firstauthor Ref. Method Organization Temperature/K Density/(g·cm-3)
    D 2000 B. Militzer [114] PIMC(VDM) University of Illinois at Urbana-Champaign 105-106 0.674,0.838
    D 2004 V. Bezkrovniy [155] DPIMC Universität Greifswald, Domstrasse, Germany;Russian Academy of Science 15 625-1 000 000 0.674,0.838,1.097
    D 2011 S. X. Hu [46] RPIMC University of Rochester,University of California 15 665-63 822 000 0.002-1 596
    H 20012003 V. S. Filinov [117]
    [116]
    DPIMC Russian Academy of Sciences 31 250-106 0.419
    H 2001 B. Militzer [114] RPIMC LLNL,University of Illinois at Urbana-Champaign 5 000-250 000 9.833×10-4-0.153
    H 2004 C. Pierleoni [121] CEIMC Universita' of L'Aquila, Via Vetoio, University of Illinois at Urbana-Champaign, Université Pierre et Marie Curie 300-10 000 5.267,2.697,1.561
    H 2010 M. A. Morales [126] CEIMC University of Illinois at Urbana-Champaign, University of L'Aquila, Italy 2 000-10 000 0.724-2.329
    H 2011 L. Caillabet [148] QMD,CEIMC CEA 116 045 0.2-5
    H 2015 Q. L. Zhang [150] DPIMC IAPCM 116 045 0.98×10-3-1 346.1
    下载: 导出CSV
  • [1] HUBBARB W B.Interiors of the giant planets[J].Science, 1981, 214(4517):145-149. doi: 10.1126/science.214.4517.145
    [2] STEVENSON D J.Interiors of the giant planets[J].Annual Review of Earth and Planetary Sciences, 1982, 10(1):257-295. doi: 10.1146/annurev.ea.10.050182.001353
    [3] MCCRORY R L, MEYERHOFER D D, BETTI R, et al.Progress in direct-drive inertial confinement fusion[J].Physics of Plasmas, 2008, 15(5):055503. doi: 10.1063/1.2837048
    [4] MEYERHOFER D D, MCCRORY R L, BETTI R, et al.High-performance inertial confinement fusion target implosions on OMEGA[J].Nuclear Fusion, 2011, 51(5):053010. doi: 10.1088/0029-5515/51/5/053010
    [5] LINDL J D.Development of the indirect-drive approach to inertial confinement fusion and the target physics basis for ignition and gain[J].Physics of Plasmas, 1995, 2(11):3933-4024. doi: 10.1063/1.871025
    [6] HU S X, GONCHAROV V N, BOEHLY T R, et al.Impact of first-principles properties of deuterium-tritium on inertial confinement fusion target designs[J].Physics of Plasmas, 2015, 22(5):056304. doi: 10.1063/1.4917477
    [7] NOH J, FULGUERAS A M, SEBASTIAN L J, et al.Estimation of thermodynamic properties of hydrogen isotopes and modeling of hydrogen isotope systems using Aspen Plus simulator[J].Journal of Industrial and Engineering Chemistry, 2017, 46:1-8. doi: 10.1016/j.jiec.2016.07.053
    [8] LEACHMAN J W, JACOBSEN R T, PENONCELLO S G, et al.Fundamental equations of state for parahydrogen, normal hydrogen, and orthohydrogen[J].Journal of Physical and Chemical Reference Data, 2009, 38(3):721-748. doi: 10.1063/1.3160306
    [9] RICHARDSON A, LEACHMAN J W, LEMMON E W.Fundamental equation of state for deuterium[J].Journal of Physical and Chemical Reference Data, 2014, 43(1):013103. doi: 10.1063/1.4864752
    [10] AZADI S, FOULKES W M C.Fate of density functional theory in the study of high-pressure solid hydrogen[J].Physical Review B, 2013, 88(1):014115. doi: 10.1103/PhysRevB.88.014115
    [11] MCMAHON J M, MORALES M A, PIERLEONI C, et al.The properties of hydrogen and helium under extreme conditions[J].Reviews of Modern Physics, 2012, 84(4):1607-1653.
    [12] KNUDSON M D, DESJARLAIS M P.High-precision shock wave measurements of deuterium:evaluation of exchange correlation functionals at the molecular-to-atomic transition[J].Physical Review Letters, 2017, 118(3):035501. doi: 10.1103/PhysRevLett.118.035501
    [13] NATOLI V V, MARTIN R M, CEPERLEY D M.Crystal structure of atomic hydrogen[J].Physical Review Letters, 1993, 70(13):1952-1955. doi: 10.1103/PhysRevLett.70.1952
    [14] SILVERA F.The solid molecular hydrogens in the condensed phase:fundamentals and static properties[J].Reviews of Modern Physics, 1980, 52(2):393-452. doi: 10.1103/RevModPhys.52.393
    [15] MAO H K, HEMLEY R J.Ultrahigh-pressure transitions in solid hydrogen[J].Reviews of Modern Physics, 1994, 66(2):671-692. doi: 10.1103/RevModPhys.66.671
    [16] MAKSIMOV E G, SHILOV Y I.Hydrogen at high pressure[J].Physics Uspekhi, 1999, 42(11):1121-1138. doi: 10.1070/PU1999v042n11ABEH000666
    [17] REDMER R, RÖPKE G.Progress in the theory of dense strongly coupled plasmas[J].Contributions to Plasma Physics, 2010, 50(10):970-985. doi: 10.1002/ctpp.201000079
    [18] GRYAZNOV V K, IOSILEVSKIY I L.Thermodynamic properties of hydrogen plasma to megabars[J].Contributions to Plasma Physics, 2016, 56(3/4):352-360.
    [19] MAO H K, CHEN B, CHEN J H, et al.Recent advances in high-pressure science and technology[J].Matter and Radiation at Extremes, 2016, 1(1):59-75. doi: 10.1016/j.mre.2016.01.005
    [20] GARBEROGLIO G, JANKOWSKI P, SZALEWICZ K, et al.Second virial coefficients of H2 and its isotopologues from a six-dimensional potential[J].The Journal of Chemical Physics, 2012, 137(15):154308. doi: 10.1063/1.4757565
    [21] SOUERS P C.Hydrogen properties for fusion energy[M].Los Angeles:University of California Press, 1986.
    [22] SAKODA N, SHINDO K, SHINZATO K, et al.Review of the thermodynamic properties of hydrogen based on existing equations of state[J].International Journal of Thermophysics, 2010, 31(2):276-296. doi: 10.1007/s10765-009-0699-7
    [23] LOUBEYRE P, LETOULLEC R, HÄUSERMANN D, et al.X-ray diffraction and equation of state of hydrogen at megabar pressures[J].Nature, 1996, 383(6602):702-704. doi: 10.1038/383702a0
    [24] HEMLEY R J, MAO H K, GONCHAROV A F, et al.Synchrotron infrared spectroscopy to 0.15 eV of H2 and D2 at megabar pressures[J].Physical Review Letters, 1996, 76:1667. doi: 10.1103/PhysRevLett.76.1667
    [25] GONCHAROV A F, GREGORYANZ E, HEMLEY R J, et al.Spectroscopic studies of the vibrational and electronic properties of solid hydrogen to 285 GPa[J].Proceedings of the National Academy of Sciences of the United States of America, 2001, 98(25):14234-14237. doi: 10.1073/pnas.201528198
    [26] DATCHI F, LOUBEYRE P, LETOULLEC R.Extended and accurate determination of the melting curves of argon, helium, ice(H2O), and hydrogen(H2)[J].Physical Review B, 2000, 61(10):6535-6546. doi: 10.1103/PhysRevB.61.6535
    [27] GREGORYANZ E, GONCHAROV A F, MATSUISHI K, et al.Raman spectroscopy of hot dense hydrogen[J].Physical Review Letters, 2003, 90(17):175701. doi: 10.1103/PhysRevLett.90.175701
    [28] NELLIS W J, MITCHELL A C, VAN THIEL M, et al.Equation-of-state data for molecular hydrogen and deuterium at shock pressures in the range 2-76 GPa(20-760 kbar)[J].The Journal of Chemical Physics, 1983, 79(3):1480-1486. doi: 10.1063/1.445938
    [29] WEIR S T, MITCHELL A C, NELLIS W J.Metallization of fluid molecular hydrogen at 140 GPa(1.4 Mbar)[J].Physical Review Letters, 1996, 76(11):1860-1863. doi: 10.1103/PhysRevLett.76.1860
    [30] KNUDSON M D, HANSON D L, BAILEY J E, et al.Principal Hugoniot, reverberating wave, and mechanical reshock measurements of liquid deuterium to 400 GPa using plate impact techniques[J].Physical Review B, 2004, 69(14):144209. doi: 10.1103/PhysRevB.69.144209
    [31] TRUNIN R F, URLIN V D, MEDVEDEV A B.Dynamic compression of hydrogen isotopes at megabar pressures[J].Physics-Uspekhi, 2010, 53(6):577-593. doi: 10.3367/UFNe.0180.201006d.0605
    [32] BELOV S I, BORISKOV G V, BYKOV A I, et al.Shock compression of solid deuterium[J].Journal of Experimental and Theoretical Physics Letters, 2002, 76(7):433-435. doi: 10.1134/1.1528696
    [33] BORISKOV G V, BYKOV A I, IL'KAEV R I, et al.Shock-wave compression of solid deuterium at a pressure of 120 GPa[J].Doklady Physics, 2003, 48(10):553-555. doi: 10.1134/1.1623535
    [34] BORISKOV G V, BYKOV A I, IL'KAEV R I.Shock compression of liquid deuterium up to 109 GPa[J].Physical Review B, 2005, 71:092104. doi: 10.1103/PhysRevB.71.092104
    [35] GRISHECHKIN S K, GRUZDEV S K, GRYAZNOV V K, et al.Experimental measurements of the compressibility, temperature, and light absorption in dense shock-compressed gaseous deuterium[J].Journal of Experimental and Theoretical Physics Letters, 2004, 80(6):398-404. doi: 10.1134/1.1830656
    [36] DA SILVA L B, CELLIERS P, COLLINS G W, et al.Absolute equation of state measurements on shocked liquid deuterium up to 200 GPa(2 Mbar)[J].Physical Review Letters, 1997, 78(3):483-486. doi: 10.1103/PhysRevLett.78.483
    [37] HICKS D G, BOEGLY T R, CELLIERS P M, et al.Laser-driven single shock compression of fluid deuterium from 45 to 220 GPa[J].Physical Review B, 2009, 79(1):014112. doi: 10.1103/PhysRevB.79.014112
    [38] SANO T, OZAKI N, SAKAIYA T, et al.Laser-shock compression and Hugoniot measurements of liquid hydrogen to 55 GPa[J].Physical Review B, 2011, 83(5):054117. doi: 10.1103/PhysRevB.83.054117
    [39] KERLEY G I.Equation of state for hydrogen and deuterium:SAND2003-3613.Sandia National Laboratories, 2003.
    [40] LOUBEYRE P, BRYGOO S, EGGERT J, et al.Extended data set for the equation of state of warm dense hydrogen isotopes[J].Physical Review B, 2012, 86(14):144115. doi: 10.1103/PhysRevB.86.144115
    [41] BORISKOV G V, BYKOV A I.Isentropic compression of substances using ultra-high magnetic field:zero isotherms of protium and deuterium the pressure range up to~5 Mbar[J].Contributions to Plasma Physics, 2011, 51(4):339-348. doi: 10.1002/ctpp.v51.4
    [42] EGOROV N I, BORISKOV G V, BYKOV A I, et al.Use of pulsed radiography for investigation of equations of state of substances at megabar pressures[J].Contributions to Plasma Physics, 2011, 51(4):333-338. doi: 10.1002/ctpp.v51.4
    [43] MOCHALOVA M A, IL'KAEVA R I, FORTOV V E, et al.Quasi isentropic compressibility of deuterium and helium at pressures of 1 500-5 000 GPa[J].Journal of Experimental and Theoretical Physics, 2014, 119(1):146-161. doi: 10.1134/S106377611406017X
    [44] FORTOV V E, ILKAEV R I, ARININ V A, et al.Phase transition in a strongly nonideal deuterium plasma generated by quasi-isentropical compression at megabar pressures[J].Physical Review Letters, 2007, 99(18):185001. doi: 10.1103/PhysRevLett.99.185001
    [45] GU Y J, CHEN Q F, ZHENG J, et al.The equation of state, shock-induced molecule dissociation, and transparency loss for multi-compressed dense gaseous H2+D2 mixtures[J].Journal of Applied Physics, 2012, 111(1):013513. doi: 10.1063/1.3675281
    [46] HU S X, MILITZER B, GONCHAROV V N, et al.First-principles equation-of-state table of deuterium for inertial confinement fusion applications[J].Physical Review B, 2011, 84(22):224109. doi: 10.1103/PhysRevB.84.224109
    [47] LIU H F, SONG H F, ZHANG Q L, et al.Validation for equation of state in wide regime:copper as prototype[J].Matter and Radiation at Extremes, 2016, 1(2):123-131. doi: 10.1016/j.mre.2016.03.002
    [48] BECKER A, LORENZEN W, FORTNEY J J, et al.Ab-initio equations of state for hydrogen(H-REOS.3)and helium(He-REOS.3)and their implications for the interior of brown dwarfs[J].The Astrophysical Journal Supplement Series, 2014, 215(2):1-21.
    [49] SAUMON D, CHABRIER G.Fluid hydrogen at high density:pressure dissociation[J].Physical Review A, 1991, 44:5122-5141. doi: 10.1103/PhysRevA.44.5122
    [50] SAUMON D, CHABRIER G.Fluid hydrogen at high density:pressure ionization[J].Physical Review A, 1992, 46:2084-2100. doi: 10.1103/PhysRevA.46.2084
    [51] LI Q, LIU H F, ZHANG G M.The thermodynamical instability induced by pressure ionization in fluid Helium[J].Physics of Plasmas, 2016, 23(11):112709. doi: 10.1063/1.4968828
    [52] MARX D, HUTTER J.Ab-initio molecular dynamics:basic theory and advanced methods[M].Cambridge, England:Cambridge University Press, 2009.
    [53] FOULKES W M C, MITAS L, NEEDS R J, et al.Quantum Monte Carlo simulations of solids[J].Reviews of Modern Physics, 2001, 73(1):33-83.
    [54] SILVERA I F, GOLDMAN V V.The isotropic intermolecular potential for H2 and D2 in the solid and gas phases[J].The Journal of Chemical Physics, 1978, 69(9):4209-4213. doi: 10.1063/1.437103
    [55] ROSS M, REE F H, YOUNG D A.The equation of state of molecular hydrogen at very high density[J].The Journal of Chemical Physics, 1983, 79(3):1487-1494. doi: 10.1063/1.445939
    [56] JURANEK H, REDMER R.Self-consistent fluid variational theory for pressure dissociation in dense hydrogen[J].The Journal of Chemical Physics, 2000, 112(8):3780-3786. doi: 10.1063/1.480939
    [57] JURANEK H, SCHWARZ V, REDMER R.Equation-of-state for hydrogen and helium in the chemical picture[J].Journal of Physics A:Mathematical and General, 2003, 36(22):6181-6185. doi: 10.1088/0305-4470/36/22/346
    [58] LIU H F, LIU W S, ZHANG W X, et al.Equations of state of H2 and D2[J].Journal of Physics Condensed Matter, 2002, 14(44):11401-11404. doi: 10.1088/0953-8984/14/44/489
    [59] ROSS M.Linear-mixing model for shock-compressed liquid deuterium[J].Physical Review B, 1998, 58(2):669-677.
    [60] CHEN Q F, CAI L C, CHEN D Q, et al.Pressure dissociation of dense hydrogen[J].Chinese Physics, 2005, 14(10):2077-2082. doi: 10.1088/1009-1963/14/10/025
    [61] LI Q, LIU H F, ZHANG G M, et al.Response to "Comment on 'The thermodynamical instability induced by pressure ionization in fluid helium'"[J].Physics of Plasmas, 2017, 24(6):064702. doi: 10.1063/1.4984999
    [62] 李琼, 刘海风, 张弓木, 等.模拟退火算法在化学自由能模型中的应用[J].计算物理, 2018.Doi:10.19596/jswl.cnki.1001-246x.2018-7853.
    [63] STOLZMANN W, BLÖCKER T.Thermodynamical properties of stellar matter Ⅰ.equation of state for stellar interiors[J].Astronomy and Astrophysics, 1996, 314:1024-1040.
    [64] BORN M, OPPENHEIMER J R.Zur Quantenthoried der Molekeln[J].Annalen der Physik, 1927, 84(20):457-484.
    [65] HARTREE D R.The wave mechanics of an atom with a non-coulomb central field.Part Ⅱ.some results and discussion[J].Mathematical Proceedings of the Cambridge Philosophical Society, 1928, 24(1):111-132. doi: 10.1017/S0305004100011920
    [66] FOCK V.Näaherungs methode zur Lösung des quanten mechanischen Mehrkörperproblems[J].Zeitschrift für Physik, 1930, 61(1/2):126-148.
    [67] Martin R M.Electronic structure:basic theory and practical methods[M].Cambridge, England:Cambridge University Press, 2004.
    [68] 谢希德, 陆栋.固体能带理论[M].上海:复旦大学出版社, 1998.
    [69] KOHN W, SHAM L J.Self-consistent equations including exchange and correlation effects[J].Physical Review, 1965, 140(4A):A1133. doi: 10.1103/PhysRev.140.A1133
    [70] KOHN W.Nobel lecture electronic structure of matter-wave functions ans density functionals[J].Reviews of Modern Physics, 1998, 71(5):1253-1266.
    [71] HERMA F, VAN DYKE J P, ORTENBURGER I B.Improved statistical exchange approxi-mation for inhomogeneous many-electron systems[J].Physical Review Letters, 1969, 22(16):807-811. doi: 10.1103/PhysRevLett.22.807
    [72] LABANOWSKI J L, ANDZELM J W.Density functional methods in chemistry[M].New York:Springer Verlag, 1991.
    [73] JUAN Y M, KAXIRAS E.Application of gradient corrections to density functional theory for atoms and solids[J].Physical Review B, 1993, 48(20):14944. doi: 10.1103/PhysRevB.48.14944
    [74] LANGRETH D C, PERDEW J P.Theory of nonuniform electronic systems.Ⅰ.analysis of the gradient approximation and a generalization that works[J].Physical Review B, 1980, 21(12):5469-5493.
    [75] PERDEW J P, BURKE K, ERNZERHOF M.Generalized gradient approximation made simple[J].Physical Review Letters, 1996, 77(18):3865-3868. doi: 10.1103/PhysRevLett.77.3865
    [76] PERDEW J P, WANG Y.Accurate and simple analytic representation of the electron gas correlation energy[J].Physical Review B, 1992, 45(23):13244. doi: 10.1103/PhysRevB.45.13244
    [77] PERDEW J P, CHEVARY J A, VOSKO S H, et al.Atoms, moleules, solids, and surfaces:applications of the generalized gradient approximation for exchange and correlation[J].Physical Review B, 1992, 46(11):6671-6687. doi: 10.1103/PhysRevB.46.6671
    [78] DION M, RYDBERG H, SCHRÖDER E, et al.Van der Waals density functional for general geometries[J].Physical Review Letters, 2004, 92(24):246401. doi: 10.1103/PhysRevLett.92.246401
    [79] LEE K, MURRAY E D, KONG L, et al.Higher-accuracy van der Waals density functional[J].Physical Review B, 2010, 82(8):081101. doi: 10.1103/PhysRevB.82.081101
    [80] IHM J, ZUNGER A, COHEN M L.Momentum-space formalism for the total energy of solids[J].Journal of Physics C:Solid State Physics, 1979, 12(21):4409-4422. doi: 10.1088/0022-3719/12/21/009
    [81] PAYNE M C, TETER M P, AHAN D C, et al.Iterative minimization techniques for ab initio total-energy calculations:molecular dynamics and conjugate gradients[J].Reviews of Modern Physics, 1992, 64(4):1045-1097. doi: 10.1103/RevModPhys.64.1045
    [82] SLATER J C, KOSTER G F.Simplified LCAO method for the periodic potential problem[J].Physical Review, 1954, 94(6):1498-1524. doi: 10.1103/PhysRev.94.1498
    [83] HERRING C, HILL A G.The theoretical constitution of metallic beryllium[J].Physical Review, 1940, 58(2):132-162. doi: 10.1103/PhysRev.58.132
    [84] HERMAN F, KILLMAN S.Atomic structure calculation[M].Englewood Cliffs, New Jersey:Prentice-Hall Inc., 1963.
    [85] ANDERSON O K.Linear methods in band theory[J].Physical Review B, 1975, 12(8):3060-3083. doi: 10.1103/PhysRevB.12.3060
    [86] SKIVER H L.The LMTO method[M].Heidelberg:Springer-Verlag, 1984.
    [87] HAMANN D R, SCHLUTER M, Chiang C.Norm-conserving pseudopotentials[J].Physical Review Letters, 1979, 43(20):1494-1497. doi: 10.1103/PhysRevLett.43.1494
    [88] BAICHELET G B, HAMANN D R, SCHLUTER M.Pseudopotentials that work:from H to Pu[J].Physical Review B, 1982, 26(8):4199-4228. doi: 10.1103/PhysRevB.26.4199
    [89] BLOCHL P E.Projector augmented wave method[J].Physical Review B, 1994, 50(24):17953-19979. doi: 10.1103/PhysRevB.50.17953
    [90] HOLZWARTH N A W, TACKETT A R, MATTHEWS G E.A projector augmented wave (PAW) code for electronic structure calculations.Part Ⅰ:atompaw for generating atom-centered functions.computer physics communications[J].Computer Physics Communications, 2001, 135(3):329-347. doi: 10.1016/S0010-4655(00)00244-7
    [91] VANDERBILT D.Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[J].Physical Review B, 1990, 41(11):7892-7895. doi: 10.1103/PhysRevB.41.7892
    [92] LAASONEN K, CAR R, LEE C, et al.Implementation of ultrasoft pseudopotentials in ab initio molecular dynamics[J].Physical Review B, 1991, 43(8):6796-6799. doi: 10.1103/PhysRevB.43.6796
    [93] KRESSE G, HAFNER J.Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements[J].Journal of Physics Condensed Matter, 1994, 6(40):8245-8257. doi: 10.1088/0953-8984/6/40/015
    [94] HELLMANN H.Einfuhrung in die Quantumchemie[M].Leipzig:Franz Duetsche, 1937.
    [95] FEYNMAN R P.Forces in molecules[J].Physical Review, 1939, 56(4):340-343. doi: 10.1103/PhysRev.56.340
    [96] MARTIN R M.Electronic structure basic theory and practical methods[J].Cambridge University Press, 2004.
    [97] CAR R, PARRINELLO M.Unified approach for molecular dynamics and density-functional theory[J].Physical Review Letters, 1985, 55(22):2471-2474. doi: 10.1103/PhysRevLett.55.2471
    [98] MARX D, HUTTE J.Modern methods and algorithms of quantum chemistry[J].NIC Series, 2000, 1:301-449.
    [99] MERMIN N D.Thermal properties of the inhomogeneous electron gas[J].Physical Review, 1965, 137(5A):A1441-A1443. doi: 10.1103/PhysRev.137.A1441
    [100] KRESSE G, HAFNER J.Ab-initio molecular dynamics for liquid metals[J].Physical Review B, 1993, 48(17):13115-13118. doi: 10.1103/PhysRevB.48.13115
    [101] GONZE X, BEUKEN J M, CARACAS R, et al.First-principles computation of material properties:the ABINIT software project[J].Computational Materials Science, 2002, 25(3):478-492. doi: 10.1016/S0927-0256(02)00325-7
    [102] GONZE X, RIGNANESE G M, VERSTRAETE M, et al.A brief introduction to the ABINIT software package[J].Zeitschrift für Kristallographie-Crystalline Materials, 2005, 220(5/6):558-562.
    [103] GONZE X, AMADON B, ANGLADE P M, et al.ABINIT:first-principles approach to material and nanosystem properties[J].Computer Physics Communications, 2009, 180(12):2582-2615. doi: 10.1016/j.cpc.2009.07.007
    [104] ABINIT. (2018-06-10). http://www.abinit.org.
    [105] GYGI F. The first-principles MD code JEEP1. 6. 6. Lawrence Livermore National Laboratory, 1999-2001.
    [106] QBOX. (2018-07-16). http://qboxcode.org.
    [107] LAMBERT F, CLÉROUIN J, MAZEVET S.Structural and dynamical properties of hot dense matter by a thomas-fermi-dirac molecular-dynamics[J].Europhysics Letters, 2006, 75(5):681-687. doi: 10.1209/epl/i2006-10184-7
    [108] MAZEVET S, LAMBERT F, BOTTIN F, et al.Ab-initio molecular-dynamics simulations of dense boron plasmas up to the semiclassical Thomas-Fermi regime[J].Physical Review E, 2007, 2007, 75:056404. doi: 10.1103/PhysRevE.75.056404
    [109] LAMBERT F, CLEROUIN J, DANEL J F, et al.Direct verification of mixing rules in the hot and dense regime[J].Physical Review E, 2008, 77:026402. doi: 10.1103/PhysRevE.77.026402
    [110] CEPERLEY D M.Path-integral calculations of normal liquid 3He[J].Physical Review Letters, 1992, 69(2):331-334. doi: 10.1103/PhysRevLett.69.331
    [111] TROTTER H.On the product of semi-groups of operators[J].Proceedings of the American Mathematical Society, 1959, 10(4):545-551. doi: 10.1090/S0002-9939-1959-0108732-6
    [112] MAGRO W R, CEPERLEY D M, PIERLEONI C, et al.Molecular dissociation in hot, dense hydrogen[J].Physical Review Letters, 1996, 76(8):1240-1243.
    [113] MILITZER B, CEPERLEY D M.Path integral monte carlo calculation of the deuterium hugoniot[J].Physical Review Letters, 2000, 85(9):1890-1893. doi: 10.1103/PhysRevLett.85.1890
    [114] MILITZER B, CEPERLEY D M.Path integral Monte Carlo simulation of the low-density hydrogen plasma[J].Physical Review, 2001, 63(2):066404.
    [115] MILITZER B, CEPERLEY D M, KRESS J D, et al.Calculation of a deuterium double shock hugoniot from ab initio simulations[J].Physical Review Letters, 2001, 87(27):275502. doi: 10.1103/PhysRevLett.87.275502
    [116] FILINOV V S, BONITZ M, LEVASHOV P R, et al.Plasma phase transition in hydrogen and electron-hole plasmas[J].Contributions to Plasma Physics, 2003, 43(5/6):290-294.
    [117] FILINOV V S, BONITZ M, EBELING W, et al.Thermodynamics of hot dense H-plasmas:path integral Monte Carlo simulations and analytical approximations[J].Plasma Physics & Controlled Fusion, 2001, 43(6):743-759.
    [118] FILINOV V S, LEVASHOV P R, BONITZ M, et al.Calculation of the shock Hugoniot of deuterium at pressures above 1 Mbar by the path-integral Monte Carlo method[J].Plasma Physics Reports, 2005, 31(8):700-704. doi: 10.1134/1.2031631
    [119] 张其黎, 刘海风, 李琼, 等.氢状态方程的路径积分蒙特卡洛研究[J].计算物理, 2018.Doi:10.19596/jswl.cnki.1001-246x.2018-7855.
    [120] PIERLEONI C, CEPERLEY D M.The coupled electron-ion Monte Carlo method[J].Lecture Notes in Physics, 2006, 703(1):641-683.
    [121] PIERLEONI C, CEPERLEY D M, HOLZMANN M.Coupled electron-ion Monte Carlo calculations of dense metallic hydrogen[J].Physical Review Letters, 2004, 93(14):146402. doi: 10.1103/PhysRevLett.93.146402
    [122] CAO J, BERNE B J.A Born-Oppenheimer approximation for path integrals with an application to electron solvation in polarizable fluids[J].The Journal of Chemical Physics, 1993, 99(4):2902-2916. doi: 10.1063/1.465198
    [123] MCMILLAN W L.Ground state of liquid He4[J].Physical Review, 1965, 138(2A):442-451. doi: 10.1103/PhysRev.138.A442
    [124] CEPERLEY D M, CHESTER G V, KALOS M H.Monte Carlo simulation of a many-fermion study[J].Physical Review B, 1977, 16(16):3081-3099.
    [125] FOULKES W M C, MITAS L, NEEDS R J, et al.Quantum Monte Carlo simulations of solids[J].Reviews of Modern Physics, 2001, 73(1):33-83. doi: 10.1103/RevModPhys.73.33
    [126] MORALES M A, PIERLEONI C, CEPERLEY D M.Equation of state of metallic hydrogen from coupled electron-ion Monte Carlo simulations[J].Physical Review E, 2010, 81(1):021202.
    [127] TUBMAN N M, LIBERATORE E, PIERLEONI C, et al.Molecular-atomic transition along the deutrium Hugoniot curve with coupled electron-ion Monte Carlo simulations[J].Physical Review Letters, 2015, 115(4):45301. doi: 10.1103/PhysRevLett.115.045301
    [128] 李名锐, 周刚, 李志康, 等.液氘单次冲击压缩的QMC模拟研究[J].高压物理学报, 2015, 29(1):1-8.
    LI M R, ZHOU G, LI Z K, et al.Single shock compression of fluid deuterium by QMC simulation[J].Chinese Journal of High Pressure Physics, 2015, 29(1):1-8.
    [129] 耿华运, 孙毅.氢的高压奇异结构与金属化[J].高压物理学报, 2018, 32(2):020101.
    GENG H Y, SUN Y.On the novel structure and metallization of hydrogen under high pressure[J].Chinese Journal of High Pressure Physics, 2018, 32(2):020101.
    [130] KERLEY G I.A theoretical equation of state for deuterium:No.LA-4776.Los Alamos, New Mexico:Los Alamos Scientific Laboratory, 1972.
    [131] SAUMON D, CHABRIER G, VAN HORN H M.An equation of state for low-mass stars and giant planets[J].The Astrophysical Journal Supplement Series, 1995, 99:713-741. doi: 10.1086/192204
    [132] EBELING W.Coulomb interaction and ionization equilibrium in partially ionized plasmas[J].Physica, 1969, 43(2):293-306. doi: 10.1016/0031-8914(69)90009-3
    [133] DA SILVA L B, CELLIERS P, COLLINS G W, et al.Absolute equation of state measurements on shocked liquid deuterium up to 200 GPa(2 Mbar)[J].Physical Review Letters, 1997, 78(78):483-486.
    [134] YOUNG D A, COREY E M.A new global equation of state model for hot dense matter[J].Journal of Applied Physics, 1995, 78(6):3748-3755. doi: 10.1063/1.359955
    [135] LENOSKY T J, KRESS J D, COLLINS L A.Molecular-dynamics modeling of the Hugoniot of shocked liquid deuterium[J].Physical Review B, 1997, 56(9):5164-5169. doi: 10.1103/PhysRevB.56.5164
    [136] GALLI G, HOOD R Q, HAZI A U, et al.Ab-initio simulations of compressed liquid deuterium[J].Physical Review B, 2000, 61(2):909-912. doi: 10.1103/PhysRevB.61.909
    [137] BAGNIER S, BLOTTIAU P, CLÉROUIN J.Local-spin-density-approximation molecular-dynamics simulations of dense deuterium[J].Physical Review E, 2001, 63(1):015301.
    [138] LENOSKY T J, BICKHAM S R, KRESS J D, et al.Density-functional calculations of the Hugoniot of shocked liquid deuterium[J].Physical Review B, 2000, 61(1):1-4.
    [139] GYGI F, GALLI G.Electronic excitations and the compressibility of deuterium[J].Physical Review B, 2002, 65(22):392-397.
    [140] DESJARLAIS M P.Density-functional calculations of the liquid deuterium Hugoniot, reshock, and reverberation timing[J].Physical Review B, 2003, 68(6):575-580.
    [141] BONEV S A, MILITZER B, GALLI G.Ab initio simulations of dense liquid deuterium:comparison with gas-gun shock-wave experiments[J].Physical Review B, 2004, 69(1):1985-1988.
    [142] DANEL J F, KAZANDJIAN L, PIRON R.Equation of state of warm dense deuterium and its isotopesfrom density-functional theory molecular dynamics[J].Physical Review E, 2016, 93(4):043210. doi: 10.1103/PhysRevE.93.043210
    [143] KARASIEV V V, CALDERÃIN L, TRICKEY S B.Importance of finite-temperature exchange correlation for warm dense matter calculations[J].Physical Review E, 2016, 93(6):063207. doi: 10.1103/PhysRevE.93.063207
    [144] KARASIEV V V, CHAKRABORTY D, SHUKRUTO O A, et al.Nonempirical generalized gradient approximation free-energy functional for orbital-free simulations[J].Physical Review B, 2013, 88(16):161108. doi: 10.1103/PhysRevB.88.161108
    [145] KNUDSON M D, DESJARLAIS M P, BECKER A, et al.Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[J].Science, 2015, 348(6242):1455. doi: 10.1126/science.aaa7471
    [146] COLLINS L A, BICKHAM S R, KRESS J D.Dynamical and optical properties of warm dense hydrogen[J].Physical Review B, 2001, 63(18):184110. doi: 10.1103/PhysRevB.63.184110
    [147] HOLST B, REDMER R, DESJARLAIS M.Thermophysical properties of warm dense hydrogen using quantum molecular dynamics simulations[J].Physical Review B, 2008, 77(18):184201. doi: 10.1103/PhysRevB.77.184201
    [148] CAILLABET L, MAZEVET S, LOUBEYRE P.Multiphase equation of state of hydrogen from ab initio calculations in the range 0.2 to 5 g/cc up to 10 eV[J].Physical Review B, 2011, 83:094101. doi: 10.1103/PhysRevB.83.094101
    [149] 刘海风, 张弓木. 液氘第一原理分子动力学模拟与实用物态方程(内部报告). 2004.
    [150] 张其黎, 张弓木, 赵艳红, 等.氘、氦及其混合物状态方程第一原理研究[J].物理学报, 2015, 64(9):094702.
    ZHANG Q L, ZHANG G M, ZHAO Y H, et al.Study of the equation of states for deuterium, helium, and their mixture[J].Acta Physica Sinica, 2015, 64(9):094702.
    [151] WANG C, ZHANG P.Wide range equation of state for fluid hydrogen from density function theory[J].Physics of Plasmas, 2013, 20(9):092703. doi: 10.1063/1.4821839
    [152] PIERLEONI C, CEPERLEY D M, BERNU B, et al.Equation of state of the hydrogen plasma by path intergral monte carlo simulation[J].Physical Review Letters, 1994, 73(16):2145-2149. doi: 10.1103/PhysRevLett.73.2145
    [153] FILINOV V S, FORTOV V E, BONITZ M, et al.Pair distribution functions of dense partially ionized hydrogen[J].Physics Letters A, 2000, 274(5):228-235.
    [154] FILINOV V S, FORTOV V E, BONITZ M, et al.Fermionic path integral Monte Carlo results for the uniform electron gas at finite temperature[J].Physical Review E, 2015, 91(3):033108. doi: 10.1103/PhysRevE.91.033108
    [155] BEZKROVNIY V, FILINOV V S, KREMP D, et al.Monte Carlo results for the hydrogen Hugoniot[J].Physical Review E, 2004, 70:057401. doi: 10.1103/PhysRevE.70.057401
    [156] RILLO G, MORALES M A, CEPERLEY D M, et al.Coupled electron-ion Monte Carlo simulation of hydrogen molecular crystals[J].Journal of Chemical Physics, 2018, 148(10):102314. doi: 10.1063/1.5001387
  • [1] 武娜田春玲刘福生匡安农袁宏宽郑兴荣 . 固氖高压物态方程的量子理论计算. 高压物理学报, 2012, 26(1): 41-47. doi: 10.11858/gywlxb.2012.01.006
    [2] 郑兴荣陈海军高晓红李继弘宋小永 . 高压固氩物态方程的量子理论计算. 高压物理学报, 2017, 31(4): 396-402. doi: 10.11858/gywlxb.2017.04.007
    [3] 龙起易范存淦邢中枢李依依 . 氢对22-13-5不锈钢室温下力学行为的影响. 高压物理学报, 1994, 8(2): 113-124 . doi: 10.11858/gywlxb.1994.02.005
    [4] 田春玲经福谦顾云军蔡灵仓刘福生 . 高温高密度氢(氘)的物态方程离解效应研究. 高压物理学报, 2007, 21(1): 8-14 . doi: 10.11858/gywlxb.2007.01.002
    [5] 张春斌 . 物态方程的分子动力学研究. 高压物理学报, 1989, 3(1): 51-57 . doi: 10.11858/gywlxb.1989.01.007
    [6] 陈其峰蔡灵仓王为傅秋卫陈栋泉经福谦 . 稠密氘气物态方程研究. 高压物理学报, 2000, 14(3): 176-181 . doi: 10.11858/gywlxb.2000.03.004
    [7] 李有宽陈栋泉李茂生董航 . 氢的高密度等离子体物态方程. 高压物理学报, 2001, 15(4): 271-276 . doi: 10.11858/gywlxb.2001.04.006
    [8] 吴绍曾张程祥吴姝妍王融朱海军 . 利用Gordon-Kim模型对NaCl物态方程的研究. 高压物理学报, 1992, 6(1): 15-22 . doi: 10.11858/gywlxb.1992.01.003
    [9] 吴强经福谦 . 用于预测疏松材料冲击压缩特性的热力学模型. 高压物理学报, 1996, 10(1): 1-5 . doi: 10.11858/gywlxb.1996.01.001
    [10] 李大红王悟 . 陶瓷物态方程实验研究. 高压物理学报, 1993, 7(3): 226-231 . doi: 10.11858/gywlxb.1993.03.010
    [11] 吴雄 . 新型爆轰产物物态方程. 高压物理学报, 1991, 5(2): 98-103 . doi: 10.11858/gywlxb.1991.02.003
    [12] 李银成 . 爆轰产物物态方程(Ⅰ)爆热、爆轰产物的等熵方程和物态方程. 高压物理学报, 1998, 12(4): 271-281 . doi: 10.11858/gywlxb.1998.04.006
    [13] 李茂生张春斌 . 关于材料稀疏区的零温物态方程. 高压物理学报, 1990, 4(1): 24-28 . doi: 10.11858/gywlxb.1990.01.004
    [14] 谢元南 . -SiC的电子结构和物态方程. 高压物理学报, 1994, 8(1): 57-59 . doi: 10.11858/gywlxb.1994.01.009
    [15] 赵艳红段素青刘海风张弓木 . RDX炸药爆轰产物物态方程. 高压物理学报, 2015, 29(1): 47-51. doi: 10.11858/gywlxb.2015.01.008
    [16] 江厚满张若棋张寿齐 . 用遗传算法确定材料物态方程参数. 高压物理学报, 1998, 12(1): 47-53 . doi: 10.11858/gywlxb.1998.01.008
    [17] 吴昆裕苏昉谢斌 . 立方晶体超声物态方程弹性参数的高压测量. 高压物理学报, 1994, 8(4): 279-284 . doi: 10.11858/gywlxb.1994.04.006
    [18] 陈其峰蔡灵仓陈栋泉经福谦赵宪庚 . 高密度流体He+H2混合物物态方程研究. 高压物理学报, 2003, 17(3): 173-179 . doi: 10.11858/gywlxb.2003.03.003
    [19] 杨信道刘福生李剑峰经福谦 . 液态甲烷的高温密度物态方程理论研究. 高压物理学报, 1998, 12(2): 120-124 . doi: 10.11858/gywlxb.1998.02.008
    [20] 孟川民姬广富黄海军 . 固氩高压物态方程和弹性性质的密度泛函理论计算. 高压物理学报, 2005, 19(4): 353-356 . doi: 10.11858/gywlxb.2005.04.012
  • 加载中
图(10)表(2)
计量
  • 文章访问数:  394
  • 阅读全文浏览量:  130
  • PDF下载量:  30
出版历程
  • 收稿日期:  2018-06-21
  • 录用日期:  2018-07-19
  • 刊出日期:  2018-10-25

氢氘物态方程研究进展

    作者简介:刘海风(1968—), 女, 硕士, 研究员, 长期从事物态方程理论研究. E-mail: liu_haifeng@iapcm.ac.cn
  • 北京应用物理与计算数学研究所, 北京 100094
基金项目:  国家自然科学基金青年项目 11604018中国工程物理研究院科学技术基金 990105国家自然科学基金重点项目 10135010科学挑战计划 TZ2016001

摘要: 针对近二十多年的氢氘物态方程理论研究工作进行了综述分析,结合本课题组改进的自由能模型、直接量子蒙特卡洛和量子分子动力学方法的模拟结果,对多个研究小组采用不同方法获得的氢氘宽区物态方程数据进行了定量评估分析。结果表明:在当前理论框架下,仅基于第一原理数值模拟得到的氢氘物态方程能够描述的热力学相空间有限;多模型集成的H-REOS.3数据库在105 K以下温度与数值模拟结果的相对差别较大,且数据稀少,二者均不能满足工程应用需求。建议采用基于半经验模型的宽区域物态方程研究方法,即结合高精度的实验研究、数值模拟和解析模型,构建满足工程应用需求的氢氘宽区实用物态方程。

English Abstract

  • 氢是自然界最丰富的元素,它及其同位素的物态方程在天体物理[1-2]、惯性约束聚变(Inertial Confinement Fusion,ICF)[3-6]、国际热核实验堆(International Thermonuclear Experimental Reactor,ITER)[7-9]等工程物理设计中具有非常重要的应用价值。氢又是元素周期表中最简单的元素,只包含一个核外电子,适于发展和检验现代凝聚态理论和数值模拟方法[10-13]。一个多世纪以来,科学家们利用多种理论和实验方法,对氢氘的物态方程及其他相关性质进行了持续的研究[14-20],特别是近二十年来,关于物态方程的半经验解析理论、先进数值模拟和实验研究工作都取得了明显进展[11-12]

    实验方面,包括静态、动态和动静耦合加载3类。早期的静态实验基于活塞圆筒等装置,主要研究压强、体积和温度(p-V-T)关系、声速及熔化特性,压力较低[21-22]。1996年以来,毛河光课题组在欧洲X射线同步辐射装置(European Synchrotron Radiation Facility,ESRF)上开展了氢氘的静态加载衍射实验,测量了固体氢在室温压缩至120 GPa的等温压缩线,还研究了固体氢压缩至285 GPa的振动谱和电子性质[23-25]。Datchi等[26]、Gregoryanz等[27]测量了氢的高压熔化曲线。动态实验中,利用气炮[28-29]Z箍缩[30]、化爆球面汇聚冲击[31-35]和高强度激光[36-38]加载,在实验样品中产生强冲击,通过测量冲击波速度、波后粒子速度,利用冲击波间断面的质量、动量和能量守恒方程,可得到雨贡纽状态的压强和密度。通常得到初始状态为常态的待测样品的冲击雨贡纽压强及密度数据,也可以通过预先压缩或加热样品、利用阻抗更高或更低的标准砧样品等方法测量更宽热力学区域的状态数据。Nellis等[28]利用气炮加载技术获得了液氘(至21 GPa)、液氢(至10 GPa)的一次冲击压缩数据,以及液氘二次冲击压缩(至76 GPa)数据。Knudson等[12, 30]利用磁驱动飞片测量了液氘21~176 GPa压力范围内的冲击雨贡纽数据,分别给出了利用Al和石英做标准材料的结果。Boriskov等[32-34]、Grishechkin等[35]利用球面汇聚压缩测量了固体、液体和气体氘60~121 GPa的一次冲击雨贡纽数据。激光冲击压缩方面:Da Silva等[36]利用Nova激光加载实验给出了液氘在25 GPa的1个和在70~210 GPa的5个冲击压缩点;Hicks等[37]利用激光加载给出了液氘45~220 GPa的一次冲击压缩数据,结果表明, 压力低于100 GPa时一次冲击压缩比小于4.2,压力在100~220 GPa之间时极限压缩比最大达5.0;Sano等[38]测量了液氢冲击加载至55 GPa的雨贡纽数据并获得了温度数据,结果与SESAME理论模型[39]一致。除了主雨贡纽外,动态冲击压缩实验还可以联合静态加压或加温技术,改变样品初始状态, 获得不同于通常初始状态的冲击压缩特性,或对样品施加多次冲击加载或柱面汇聚加载, 获得准等熵压缩数据。Loubeyre等[40]利用动-静态联合实验,即利用金刚石压砧(Diamond Anvil Cell,DAC)装置预先压缩氢样品,然后用激光加载测量氢的冲击压缩特性。Weir等[29]利用样品和窗口间的冲击波多次反射,测量了液氢的多次冲击压缩数据。Boriskov等[41-44]利用爆磁压缩装置测量了氢氘的准等熵数据。国内陈其峰等[45]利用二级轻气炮加载方式测量了氢氘混合物的状态方程等。可信的状态方程数据测量对冲击波速度的测量精度提出了很高的要求,同时数据分析也需要基于多种假设,如样品初始状态为常态或状态量能被准确测量、冲击波后物质处在平衡态、样品仅受到一维压缩且沿径向密度均匀等,相比较而言,对于以上假设主雨贡纽的测量更易满足,偏离雨贡纽实验分析需要考虑更复杂的实验因素。

    上述实验研究提供了较为丰富的状态方程数据。但是,由于ICF等工程物理过程中,氘氚经历的热力学范围非常广,涉及的物质状态包括固、液、气及等离子体状态[46-48],在不同的状态下,物质基本组成成分复杂,包括分子、原子、分子离子、原子离子及电子。在温度从临界温度到上亿摄氏度、密度从临界密度到几千克每立方厘米的热力学相空间内,离子的关联相互作用有强有弱,电子的存在形式包括从束缚到自由及众多中间状态。实际工程可用的宽区理论物态方程研究极其困难,需要结合实验数据,通过一定的理论框架进行合理延拓,对于不具备解析理论基础或解析理论准确度不够的热力学区域,也发展了定量的数值模拟方法,得到了大量的理论数据,但这些理论物态方程的准确度如何仍没有明确的判断,造成了工程设计中关于物态方程选取的困惑,也成为写作本文的主要动力。

    本文针对近二十多年的氢氘物态方程理论研究工作进行了综述分析,结合课题组改进的自由能模型、直接量子蒙特卡洛(Direct Path-Integral Monte Carlo,DPIMC)和量子分子动力学(First-Principles Molecular Dynamics,FPMD; 或Quantum Molecular Dynamics,QMD)方法的模拟结果,对不同方法获得的大量物态方程数据进行定量评估分析。最后建议采用多模型优化集成的技术路线,结合高精度的实验研究、数值模拟和解析模型,构建满足工程应用要求的氢氘宽区实用物态方程。

    • ICF靶丸压缩过程中,氘氚经历的温度、密度范围覆盖了从低温低密度到高温高密度的宽广热力学区域,不同的温度、密度区域对应不同的物理问题,因而需要采用不同的物理模型和理论方法。通常而言,理想气体模型适用于密度足够低且温度适中的区域;当温度降低时,气体分子的量子性逐渐表现出来,此区域要用量子理想气体模型,其物态方程也具有简明的解析形式。当密度不太高、温度升高时,会发生分子离解和原子电离,常称此为温致离化,该区域的粒子组分和物态方程要用萨哈(Saha)理论计算。当温度不太高、密度增大时,需要考虑粒子之间的相互作用,如:对于由分子组成的气体和液体,粒子之间为短程相互作用,可用某种函数描述相互作用势,然后采用多种形式的级数展开(如分子集团积分展开方法)计算物态方程,也可以直接采用经验物态方程(如范德瓦尔斯方程);对于稠密气体和液体,常用液体变分微扰理论。当温度不太高、密度增大时,也会发生分子离解和原子电离,即所谓压致离化。在温度或密度极高的区域,温致离化或压致离化彻底完成,称为完全离化区。完全离化区的低密度部分可用理想气体模型,高密度部分可用托马斯费米(Thomas-Fermi,TF)及其修正模型,也可用以量子多体理论和蒙特卡洛模拟为基础的完全电离等离子体模型。对于某些温度和密度区域,通常同时存在温致离化和压致离化,系统成分包括分子、原子、各级离子和电子,常用的模型可分为两类:一类基于化学图景(the chemical picture),视物质的基本构成粒子为分子、原子、离子和公有化电子;另一类基于物理图景(the physical picture),视物质的基本构成粒子为电子和原子核。基于化学图景的模型与传统的Saha、液体微扰论类似,在其原有基础上,增加系统的离子组分种类,扩展自由能的表述方式,构成广义的化学模型或自由能模型,本文通称此类模型为自由能模型[49-51]。基于物理图景的模型从求解电子和核构成的多体系统哈密顿方程出发,电子和离子之间通过库仑势相互作用,用量子理论描述原子和分子的形成,通过大规模的数值模拟获得系统的性质,常用的方法包括QMD方法[52]和量子蒙特卡洛(Quantum Monte Carlo,QMC)方法[53]

      关于极端条件下氢和氦的性质研究情况,McMahon等[11]作了详细综述。为便于理解后面的物态方程定量分析结果,本文从自由能模型、量子分子动力学模拟、量子蒙特卡洛模拟、实用宽区域理论模型4个方面进行概要介绍。

    • 当温度、密度足够高时,原子完全离化,等离子体模型和屏蔽库仑势模型可以较好地描述等离子体性质。对部分离化区域,可假定系统由一组确定的化学组分构成,组分之间的相互作用近似表示为某种解析形式,其参数可以利用某些特定实验数据或第一原理结果确定;同时假定系统配分函数包括内部(振动转动和电子激发)和外部(质心相互作用)自由度,直接写出系统的自由能,通过求解自由能极小值确定系统组分,并计算系统的其他热力学量。Silvera等[54]最早用这类模型计算低密度氢的物态方程,系统仅含氢分子,氢分子相互作用采用SG(Silvera-Goldman)势。为使这个模型可用于更宽的温度和压力区域,Ross等[55]、Juranek等[56-57]开展了更细致的组分或相互作用势等的改进工作,给出了适用于不同区域的物态方程数据:Ross等[55]采用软化短程相互作用的势函数,复现了更高压力下的实验数据;Juranek等[56]采用液体变分微扰理论研究了氢分子的离解平衡;Juranek等[57]采用化学模型研究了氢分子的离解和电离。国内刘海风等[58]基于理想混合模型[59]计算了氢氘的低压雨贡纽和等温线;陈其峰等[60]利用自洽变分液体微扰模型,研究了氢在密度0.01~1.0 g/cm3、温度2 000~10 000 K范围内的离解度和物态方程。

      最近,李琼等[51, 61-62]根据化学模型思想建立了改进的自由能模型,对氢的宽区域物态方程进行了系统研究,详细模型将另文总结,这里仅介绍梗概。该模型中,氢在部分离化区被视为H2、H、H+和e的混合物,粒子组分由热力学第二定律——自由能函数取极小值决定。混合物的自由能包含各种粒子的平动能、相互作用能和内部运动能量。相互作用能进一步分为3部分:中性粒子(H2、H)之间的范德瓦尔斯相互作用能、带电粒子(H+、e)之间的库仑相互作用能、中性粒子与带电粒子之间的电偶极化相互作用能。平动能部分按如下方式计算:“重粒子”H2、H、H+遵循麦克斯韦-玻耳兹曼统计,电子e遵循费米-狄拉克统计。中性粒子(H2、H)之间的范德瓦尔斯相互作用势用EXP-6势函数描述。对于多组元混合物,可采用混合硬球变分液体微扰理论处理,也可基于等效单组元近似采用单纯硬球变分液体微扰理论处理。等效单组元近似把多组元混合物视为一种等效组元的纯净物,等效组元的相互作用势与混合组分相关。库仑相互作用能采用完全电离等离子体模型的Pade近似多项式[63]计算。带电粒子的库仑势导致中性粒子发生电偶极化,采用维里展开计算电偶极化相互作用能[50]。内部运动能量包含氢分子的振动转动能和氢原子内部的电子运动能量。李琼考虑了振动的非谐效应、转动的形变效应以及振动和转动的耦合效应,应用SESAME[39]提供的数值拟合公式计算分子的转动振动自由能。此外,还分析了相互作用势、电子能量截断等多种因素对物态方程数据的影响。改进模型对体系的自由能作了更细致、准确的描述,大大拓宽了自由能模型适用的温度、密度范围。并且自由能模型需要的计算量相对较小,可清楚地给出各成分对物态方程的贡献,在宽区域物态方程的构建中将发挥重要作用。

    • 实际材料体系的主要物理性质本质上都取决于电子结构,而电子结构由其波函数确定。薛定谔方程提供了对多电子系统微观粒子运动的统一描述

      $ \mathit{\boldsymbol{H\psi }}\left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{R}}} \right) = E\mathit{\boldsymbol{\psi }}\left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{R}}} \right) $

      式中:Ψ(r, R)为体系的波函数,rR分别代表所有电子和原子核的坐标,H为体系的哈密顿量。通过求解薛定谔方程,获得系统的波函数,就可以得到系统的其他物理量。体系的哈密顿量由3部分组成,分别是电子能量He、原子核能量HN、电子与核的相互作用能He-N,即

      $ \begin{array}{l} \mathit{\boldsymbol{H}} = {\mathit{\boldsymbol{H}}_{\rm{e}}} + {\mathit{\boldsymbol{H}}_{\rm{N}}} + {\mathit{\boldsymbol{H}}_{{\rm{e - N}}}} = \\ \;\;\;\;\;\;\; - \sum\limits_i {\frac{{{\hbar ^2}}}{{2{m_{\rm{e}}}}}\nabla _{{r_i}}^2} + \frac{1}{2}\sum\limits_{i \ne i'} {\frac{{{e^2}}}{{\left| {{\mathit{\boldsymbol{r}}_i} - {\mathit{\boldsymbol{r}}_{i'}}} \right|}}} - \sum\limits_j {\frac{{{\hbar ^2}}}{{2{M_j}}}\nabla _{{R_j}}^2} + \frac{1}{2}\sum\limits_{j \ne j'} {{V_{\rm{N}}}\left( {{\mathit{\boldsymbol{R}}_j} - {\mathit{\boldsymbol{R}}_{j'}}} \right)} - \sum\limits_{i,j} {{V_{{\rm{e}} - {\rm{N}}}}\left( {{\mathit{\boldsymbol{r}}_j} - {\mathit{\boldsymbol{R}}_j}} \right)} \end{array} $

      真实材料是包含阿伏伽德罗常数量级电子和原子核的多粒子系统,其系统哈密顿量非常复杂,无法直接求解薛定谔方程,必须采取一些近似和简化处理,首先通过波恩-奥本海默绝热近似(Born-Oppenheimer approximation,BO)[64]将原子核和电子的运动分开,将研究对象抽象为具有平移周期对称的理想体系,再通过哈特利-福克自洽场[65-66]方法将多电子问题简化为单电子问题,然后通过密度泛函理论[67-69](Density Functional Theory,DFT)实现求解。

      单电子基态自洽Kohn-Sham(KS)方程[69-70]形式上类似于Schrödinger方程,包含了所有的多体效应

      $ {E_{{\rm{KS}}}}\left[ n \right] = - \frac{1}{2}\sum\limits_{i = 1}^N {{{\left| {\vec \nabla {\psi _i}\left( \mathit{\boldsymbol{r}} \right)} \right|}^2}} + \int {{{\rm{d}}^3}\mathit{\boldsymbol{r}}n\left( \mathit{\boldsymbol{r}} \right){V_{{\rm{ext}}}}\left( \mathit{\boldsymbol{r}} \right)} + {E_{\rm{H}}}\left[ n \right] + {E_{{\rm{NN}}}} + {E_{{\rm{xc}}}}\left[ n \right] $

      $ {E_{\rm{H}}}\left[ n \right] = \frac{1}{2}\int {\int {{{\rm{d}}^3}\mathit{\boldsymbol{r}}{{\rm{d}}^3}\mathit{\boldsymbol{r'}}\frac{{n\left( \mathit{\boldsymbol{r}} \right)n\left( {\mathit{\boldsymbol{r'}}} \right)}}{{\left| {\mathit{\boldsymbol{r}} - \mathit{\boldsymbol{r'}}} \right|}}} } $

      式中:$ n\left( \mathit{\boldsymbol{r}} \right) = \sum\limits_{i = 1}^N {{{\left| {{\psi _i}\left( \mathit{\boldsymbol{r}} \right)} \right|}^2}} $是无相互作用系统电子密度,ψi(r)是无相互作用系统Hamiltonian的本征态, EH是Hatree能,ENN是核-核相互作用能,Exc[n]是交换关联能。在KS形式中,所有多体效应都包含在交换关联泛函Exc[n]中,将求解基态电荷密度的多电子系统问题在形式上用描述单电子运动的等效KS方程组代替,即相互作用多体系统的基态问题在形式上严格地转化为有效势场中运动的独立粒子基态问题,有效势场由系统中所有电子贡献自洽决定。

      在KS方程中,交换关联泛函Exc[n]没有明确的形式,近似主要集中在交换关联项,交换关联泛函的准确性决定了密度泛函理论的精度。最常见、最简单有效的近似是局域密度近似(Local Density Approximation,LDA)[71-73]及在其基础上改进的广义梯度近似(Generalized Gradient Approximation,GGA)[74-75]。GGA包含密度梯度展开的能量泛函,有很多不同的表述形式,目前使用最广泛的GGA泛函有Perdew-Burke-Ernzerhof(PBE)[75]形式,以及Perdew-Wang(PW91)[76-77]、PBE96[75]等。

      基于LDA、GGA近似的密度泛函理论构成了目前最流行的电子结构计算方案,但针对氢氘分子、原子的转变问题,改进的泛函(如非局域交换关联函数vdW-DF1[78]、vdW-DF2[79])的模拟给出了与实验更一致的结果。

      交换-关联泛函形式确定后就完全定义了能量泛函,KS方程给出了与Hartree-Fock方程相似的单电子方程,求解该方程要选择一组合适的完备集基函数将本征函数展开,即把KS方程表示为矩阵形式,求矩阵的本征值和本征矢。选取合适的基函数在自洽过程中非常重要。原则上基函数应包含完备基中所有的基矢,但当基函数个数很大时,久期方程维度太高,计算量很大,所以需尽可能包含少的基函数。根据不同的研究对象,有很多基函数选择方法,如赝势平面波(Pseudopotential Plane-Wave)法[80-81]、原子轨道线性组合(Linear Combination of Atomic Orbitals,LCAO)[82]、正交平面波(Orthogonalized Plane Wave,OPW)法[83]、缀加平面波(Augmented Plane Wave,APW)法[84]、线性缀加平面波(Linearized Augmented Plane Wave,LAPW)法[85]、Mufin-Tin轨道线性组合方法(Linear Muffin-Tin Orbital Method,LMTO)[86]等。

      对周期边界条件,平面波基矢是很自然的选择,它表示扩展的、非局域态, 是最简单的正交完备基。根据Bloch定理,单电子波函数可以用平面波展开

      $ {\psi _i}\left( \mathit{\boldsymbol{r}} \right) = \sum\limits_\mathit{\boldsymbol{G}} {{c_i}\left( \mathit{\boldsymbol{G}} \right){{\rm{e}}^{{\rm{i}}\mathit{\boldsymbol{G}} \cdot \mathit{\boldsymbol{r}}}}} $

      式中:G是原胞的倒格矢,求和遍及所有|G|2/2<Ecutoff的倒格矢;Ecutoff是截断能,原则上构成完备集需要无穷多个平面波,但具有较小动能|K+G|2/2的平面波系数ci(K+G)比具有较大动能的平面波系数大,因此可以只用小于某一能量Ecutoff的平面波作为基进行展开。Ecutoff越小,计算量越小,但截断所引起的误差越大,因而需要增加Ecutoff直到收敛,Ecutoff的选择取决于赝势以及电子态的性质。交换-关联能很容易在实空间网格中计算,而动能和库仑能在倒空间求和,在求解过程中,需要通过快速傅里叶变换(Fast Fourier Transform,FFT)在实空间和倒空间计算。

      固体能带计算实际求解KS方程时,在原子核附近,芯态电子的轨道有很强的局域性,振荡剧烈,如采用平面波作为基矢来描述,则需要非常大的基组,将给计算带来很大的困难。在绝大多数条件下,芯态电子不参与成键,固体的热力学和化学性质主要由费米能级附近的价态电子决定。价电子在材料的特性上扮演极重要的角色,参与了电荷转移及成键,基于此,在计算中将芯态电子和原子核当作离子实整体考虑,其对价电子的作用由一个等效势描述,通常将该有效势称为赝势,利用赝势得到的本征能量与真实本征值一致,波函数在离子实外与真实波函数一致,在离子实内缓慢变化。赝势包含了核的库仑吸引和芯电子的屏蔽效应,是重整化的关于价态的电子-核相互作用势。在同样的计算精度下,赝势的总能计算减少了平面波基矢数目,效率比全电子方法高得多。

      赝势的导出不是唯一的,目前有很多种类的赝势,它们具有不同水平的可移植性及复杂性。最常用的赝势是Hamann提出的模守恒赝势[87-88]、Vanderbilt提出的超软赝势[89-91]以及投影缀加波(Projector Augmented Wave,PAW)[92-93]方法。模守恒赝势所对应的波函数在芯区以外(rrc)和真实波函数的形状及幅度相同,所谓模守恒是通过保证芯区(rrc)内赝波函数和真实波函数模平方求积分一致而实现的,模守恒赝势可用于不同的化学环境,移植性较好,但由于定域性较强,需要平面波能量截断很高。Vanderbilt提出的超软赝势[89-91]放弃了模守恒条件,芯区赝波函数十分光滑,根据广义的正交条件,可以采用尽可能小的平面波截断能,Vanderbilt超软赝势降低了计算量,在实际中应用广泛。

      PAW方法[92-93]包含了赝势和全电子方法的优点,计算量却没有明显增加,其基本思想是把以原子为中心的缀加球分成两部分,分别对待离原子核较近和较远的波函数:离原子核远的波函数由于其变化平缓,PAW方法不做改变; 把原子核周围近距离的波函数变换到一个赝Hilbert空间中,在这个赝空间处理波函数,再转换回原空间。这样就在赝Hilbert空间建立了赝波函数的KS方法,得到基态。缀加平面波方法也有很多近似处理,如冻结芯电子、分波截断等。在实际高压物性研究应用中,产生的PAW数据要保持缀加球半径足够小,不能重叠,需要的平面波截断能也随着半径的减小而增大。

      经典分子动力学忽略了电子的运动,没有考虑量子效应,只考虑了原子核的运动。原子核的运动遵从牛顿第二定律,相互作用力一般来自经验势函数,而经验势函数的获取一般都来自实验值,当温度、压强超出实验范围较远时,经验势就变得不可靠。第一原理分子动力学从电子波函数出发,釆用Hellmann-Feynman[94-96]定理获得离子间的相互作用力,离子运动仍然遵循牛顿第二定律,理论上可实现较宽热力学区域的物性模拟。

      1985年,Car和Parrinello[97]成功地将分子动力学和DFT有机地结合起来,提出了第一原理分子动力学方法,即Car-Parrinello分子动力学(Car-Parrinello Molecular Dynamics,CPMD),扩展了第一原理计算方法模拟的领域。CPMD采用BO近似,电子KS轨道和离子的动力学同时进行,用拉格朗日(Lagrangian)函数定义动力学与势能曲面有关,而势能曲面不但是离子位置的函数,还是电子自由度的函数EKS=EKS(R, {ψi}),用Lagrangian乘子保证KS轨道的正交性,可写为

      $ L = \sum\limits_{i = 1}^{{N_{\rm{e}}}} {\mu \int {{\rm{d}}\mathit{\boldsymbol{r}}{{\left| {{{\dot \psi }_i}\left( \mathit{\boldsymbol{r}} \right)} \right|}^2}} } + \sum\limits_{I = 1}^{{N_{\rm{n}}}} {\frac{1}{2}{M_I}\mathit{\boldsymbol{\dot R}}_I^2} - {E^{{\rm{KS}}}}\left( {\psi ,\mathit{\boldsymbol{R}}} \right) + \sum\limits_{ij} {{\mathit{\Lambda }_{ij}}\left[ {\left\langle {{\psi _i}\left| {{\psi _j}} \right.} \right\rangle - {{\rm{ \mathsf{ δ} }}_{ij}}} \right]} $

      式中:MIRI$ {\mathit{\boldsymbol{\dot R}}_I} $分别为离子质量、坐标和速度,ψi是电子波函数,EKS为对应离子位置RI时系统的Kohn-Sham能量,μ是赋予电子自由度的虚质量。由此得到电子轨道和离子位置的运动方程

      $ \mu {{\ddot \psi }_i}\left( {\mathit{\boldsymbol{r}},t} \right) = - {{\hat H}_{{\rm{KS}}}}{\psi _i}\left( {\mathit{\boldsymbol{r}},t} \right) + \sum\limits_k {{\mathit{\Lambda }_{ik}}{\psi _k}\left( {\mathit{\boldsymbol{r}},t} \right)} $

      $ {M_I}\mathit{\boldsymbol{\ddot R}} = \mathit{\boldsymbol{F}}_I^{{\rm{KS}}} = - \partial {E^{{\rm{KS}}}}/\partial {\mathit{\boldsymbol{R}}_I} $

      利用使KS能量泛函极小的电子轨道开始模拟,离子位置和电子轨道在模拟中同时演化。调节虚质量使电子子系统尽量接近BO能量曲面,同时离子保持在物理温度T。要成功地应用该方法,在整个模拟过程中电子自由度都应该在BO能量曲面附近振荡。这只有在如下条件下才能做到:电子和离子自由度之间没有能量流动,或者能流足够小使得在电子的热效应变得明显前可以进行足够长时间的模拟。要满足轨道随离子绝热移动,不发生能量转移,两个子系统的功率谱必须在频率范围内没有重叠。在有电子能隙的系统,仔细调节电子轨道的虚质量可以做到这一点。

      另一种最简单、最直接的方法是波恩-奥本海默分子动力学(Born Oppenheimer Molecular Dynamics,BOMD)[98]。BOMD方法中,系统的拉格朗日函数为

      $ {\ell _{{\rm{BO}}}}\left[ {{\mathit{\boldsymbol{R}}_I},{{\mathit{\boldsymbol{\dot R}}}_I}} \right] = \sum\limits_{I = 1}^{{N_{\rm{n}}}} {\frac{1}{2}{M_I}\mathit{\boldsymbol{\dot R}}_I^2} - \mathop {\min }\limits_{\left\{ {{\psi _i}} \right\}} {E^{{\rm{KS}}}}\left[ {\left\{ {{\psi _i}} \right\},{\mathit{\boldsymbol{R}}_I}} \right] $

      离子运动方程仍为(8)式,但由于计算中考虑了温度效应,实际用于变分的系统能量不是(3)式中的EKS,而是基于Mermin函数的有限温度密度泛函自由能[99]

      $ F\left[ {n\left( r \right),{f_i},\mu } \right] = E\left[ {n\left( r \right)} \right] - \sigma s\left[ {{f_i}} \right] + \mu \left[ {\sum\limits_i {{f_i}} - {N_{\rm{e}}}} \right] $

      式中:fi是电子占据数;μ是电子化学势;σ=kBTkB是Boltzmann常数;s[fi]是无相互作用系统的电子熵;Ne为总电子数。BOMD方法在分子动力学的每一步都进行完全收敛的DFT计算,通过求解有限温度能量泛函极小时的轨道和电子密度,根据Hellman-Feynman(HF)计算离子受力。BOMD模拟要求完全收敛的DFT计算,非常耗费机时,但随着快速高效的DFT算法的发展,BOMD逐渐比CPMD使用得更广泛,因为它更容易控制,且不容易出错。

      根据这两种不同的第一原理分子动力学物理方案,多个研究团体发展了多种计算程序,如常见的VASP[100-103](维也纳科技大学开发的第一原理模拟程序)、ABINIT[104]、JEEP(CPMD方法)[105]、QBOX[106]等程序。这些程序从根本上讲都是基于第一性原理的数值模拟程序,但由于实现数值计算的方法不同,应用程序时关于计算体系或超原胞大小的选择、K点取样、电子波函数的平面波基矢截断能量、计算时间、交换关联函数等选择不同,可能导致计算的物态方程有明显差别。

      由于计算机硬件的限制,基于量子分子动力学的数值模拟工作在2000年后才得到深入开展,此方法用于氢及氘的物态方程模拟工作非常多,详见2.1节。

      在电子结构的自洽计算中,QMD方法需要求解与单粒子轨道相关的KS方程,随着温度的升高,需要计算的轨道数目迅速增加,导致计算量大大增加,无轨道分子动力学方法(Orbits Free Molecular Dynamics, OFMD)[107-109]可用于高温区的计算。但由于OFMD方法中动能项通过TF或其修正模型表示为密度的泛函,物理实质是TF及其修正模型的数值化,同时考虑了离子的关联,因此模拟精度与TF等模型处在同一水平,后面的数据评估中不再关注该方法的模拟结果。

    • 量子多体系统在有限温度下的物理量期望值可表示为

      $ {\left\langle {\hat O} \right\rangle _\rho } = {Z^{ - 1}}\left\langle {\hat O\hat \rho } \right\rangle = {Z^{ - 1}}\int {{\rm{d}}\mathit{\boldsymbol{R}}{\rm{d}}\mathit{\boldsymbol{r}}{\rm{d}}\mathit{\boldsymbol{R'}}{\rm{d}}\mathit{\boldsymbol{r'}}\rho \left( {\mathit{\boldsymbol{R}},\mathit{\boldsymbol{r}};\mathit{\boldsymbol{R'}},\mathit{\boldsymbol{r'}};\beta } \right)\left\langle {\mathit{\boldsymbol{R'}},\mathit{\boldsymbol{r'}}\left\langle {\hat O} \right\rangle \mathit{\boldsymbol{R}},\mathit{\boldsymbol{r}}} \right\rangle } $

      式中:β=1kBTT是温度,Z是配分函数

      $ Z = \int {{\rm{d}}\mathit{\boldsymbol{R}}{\rm{d}}\mathit{\boldsymbol{r}}\rho \left( {\mathit{\boldsymbol{R}},\mathit{\boldsymbol{r}};\mathit{\boldsymbol{R'}},\mathit{\boldsymbol{r'}};\beta } \right)} $

      ρ(R, r; R', r'; β)是给定系综下在实空间表象的密度矩阵

      $ \rho \left( {\mathit{\boldsymbol{R}},\mathit{\boldsymbol{r}};\mathit{\boldsymbol{R'}},\mathit{\boldsymbol{r'}};\beta } \right) = \left\langle {\mathit{\boldsymbol{R}},\mathit{\boldsymbol{r}}\left| {{{\rm{e}}^{ - \beta \hat H}}} \right|\mathit{\boldsymbol{R'}},\mathit{\boldsymbol{r'}}} \right\rangle $

      根据exp[-(β1+β2)$\hat H$]=exp(-β1$\hat H$)exp(-β2$\hat H$),可将(12)式写成路径积分的形式

      $ Z = \int {\prod\limits_{t = 0}^{M - 1} {{\rm{d}}{\mathit{\boldsymbol{R}}^t}{\rm{d}}{\mathit{\boldsymbol{r}}^t}\left\langle {{\mathit{\boldsymbol{R}}^t},{\mathit{\boldsymbol{r}}^t}\left| {{{\rm{e}}^{ - \tau \hat H}}} \right|{\mathit{\boldsymbol{R}}^{t + 1}},{\mathit{\boldsymbol{r}}^{t + 1}}} \right\rangle } } $

      式中:τ=β/M。对可分辨粒子,沿t方向应满足周期边界条件:R0=RMr0=rM。对满足Bose或Fermi统计的不可分辨粒子,还要考虑粒子之间的交换对称性。Ceperley[110]利用Trotter分解公式[111]将密度矩阵表示为如下形式

      $ \rho \left( {{\mathit{\boldsymbol{R}}^0},{\mathit{\boldsymbol{r}}^0};{\mathit{\boldsymbol{R}}^0},{\mathit{\boldsymbol{r}}^0};\beta } \right) = \frac{1}{{N!}}\sum\limits_{\rm{P}} {{{\left( { \pm 1} \right)}^{\rm{P}}}\int_{{\mathit{\boldsymbol{R}}^0} \to {\rm{P}}{\mathit{\boldsymbol{R}}^0}} {{\rm{d}}\mathit{\boldsymbol{X}}{{\rm{e}}^{ - S\left( \mathit{\boldsymbol{X}} \right)}}} } $

      式中:S表示路径的作用量,X={R0, r0, …, RM, rM}代表路径,包含3M(Ne+Nn)个变量,NeNn分别表示电子和原子核数;P表示粒子之间的交换,公式中对所有的交换求和,对满足Bose统计的交换取“+”号,满足Fermi统计的交换取“-”号。

      对电子和原子核组成的系统,配分函数形式可表示为[11]

      $ Z = \oint {{\rm{D}}\mathit{\boldsymbol{R}}{\rm{exp}}\left[ { - \int_0^\beta {{\rm{d}}t{T_{\rm{n}}}\left( {{\mathit{\boldsymbol{R}}^t}} \right)} } \right]{Z_{{\rm{el}}}}\left( {{\mathit{\boldsymbol{X}}_{\rm{n}}}\left( t \right),\left[ {{\mathit{\boldsymbol{X}}_{\rm{n}}}} \right]} \right)} $

      $ {Z_{{\rm{el}}}}\left( {{\mathit{\boldsymbol{X}}_{\rm{n}}}\left( t \right),\left[ {{\mathit{\boldsymbol{X}}_{\rm{n}}}} \right]} \right) = \oint {{\rm{D}}\mathit{\boldsymbol{r}}{\rm{exp}}\left[ { - \int_0^\beta {{\rm{d}}t{S_{{\rm{el}}}}\left( {{\mathit{\boldsymbol{R}}^t},{\mathit{\boldsymbol{r}}^t}} \right)} } \right]} $

      式中:∮DR(∮Dr)是对所有路径求泛函积分的缩写,Xn表示路径X的原子核坐标,Tn表示原子核的动能,Zel是路径中原子核坐标的函数。

      对Bose子(取“+”号)或Boltzmann粒子(没有粒子之间的交换求和),(15)式是正定的,其配分函数与由(Ne+Nn)M个满足Boltzmann分布的粒子组成的经典系统相同,粒子之间通过有效的经典“势”kBTS(X)相互作用。对电子系统,由于交换反对称性,(15)式有正有负,配分函数的积分计算由于正负相消而变得复杂,这就是所谓的“负符号”问题。实际路径积分计算方案中,衍生了多种计算方法。常用的方法有3类[11]:第1类是非零温系统的路径积分蒙特卡洛方法,包括约束路径积分(Restricted Path-Integral Monte Carlo,RPIMC)和直接费米路径积分蒙特卡洛(Direct Path-Integral Monte Carlo,DPIMC)方法;第2类采用波恩奥本海默近似,将电子和核的运动分离,采用MC方法处理有限温度核的运动,采用零温基态QMC方法处理电子运动,发展了耦合电子-离子蒙特卡洛(Coupled Electron-Ion Monte Carlo,CEIMC)方法;第3类主要针对零温费米子系统,波函数为基的方法包括变分蒙特卡洛(Variational Monte Carlo,VMC)、爬虫蒙特卡洛(Reptation Monte Carlo, RMC)、扩散蒙特卡洛(Diffusion Monte Carlo,DMC)方法。

      为了解决“负符号”问题,Ceperley[110]发展了RPIMC方法:在路径上引入参考点R*,以确定密度矩阵相对于R*的节点ρ(R, R*, t)=0,在计算中只对避开节点的路径进行积分,即:在0<tβ, ρ(R(t), R*, t)≠0。通过约束

      $ {\rho _{\rm{F}}}\left( {{\mathit{\boldsymbol{R}}_\beta },{\mathit{\boldsymbol{R}}^ * },\beta } \right) = \int {{\rm{d}}{\mathit{\boldsymbol{R}}_0}{\rho _{\rm{F}}}\left( {{\mathit{\boldsymbol{R}}_0},{\mathit{\boldsymbol{R}}^ * },0} \right)} \oint_{{\mathit{\boldsymbol{R}}_0} \to {\mathit{\boldsymbol{R}}_\beta } \in Y\left( {{\mathit{\boldsymbol{R}}^ * }} \right)} {{\rm{d}}{\mathit{\boldsymbol{R}}_t}{{\rm{e}}^{ - S\left[ {{\mathit{\boldsymbol{R}}_t}} \right]}}} $

      所有积分的路径对配分函数的贡献都是正值。原则上,如果在约束条件中用到的密度矩阵是准确的,则符号问题能得到真正的解决; 但实际上,准确的密度矩阵是未知的,通常只能采用某种近似。最简单的近似是采用单粒子密度矩阵的行列式

      $ \rho \left( {\mathit{\boldsymbol{R}},\mathit{\boldsymbol{R'}};\beta } \right) = \left| {\begin{array}{*{20}{c}} {{\rho _1}\left( {{\mathit{\boldsymbol{r}}_1},{{\mathit{\boldsymbol{r'}}}_1};\beta } \right)}& \cdots &{{\rho _1}\left( {{\mathit{\boldsymbol{r}}_N},{{\mathit{\boldsymbol{r'}}}_1};\beta } \right)}\\ \vdots&\vdots&\vdots \\ {{\rho _1}\left( {{\mathit{\boldsymbol{r}}_1},{{\mathit{\boldsymbol{r'}}}_N};\beta } \right)}& \cdots &{{\rho _1}\left( {{\mathit{\boldsymbol{r}}_N},{{\mathit{\boldsymbol{r'}}}_N};\beta } \right)} \end{array}} \right| $

      单粒子密度矩阵可采用自由粒子密度矩阵或超越自由粒子的密度矩阵。要注意的是,在这些方法中,试探密度矩阵仅用于对节点的约束,而在路径积分作用量的计算中则要考虑完整的相互作用势。RPIMC方法已经用于对氢(氘)的物性研究[11, 112-115]

      Filinov等[116-118]发展了DPIMC方法。对于由Ne个电子和Nn个质子组成的二元混合系统,其配分函数Z可表示为如下形式

      $ Z\left( {{N_{\rm{e}}},{N_{\rm{n}}},V,\beta } \right) = \frac{{Q\left( {{N_{\rm{e}}},{N_{\rm{n}}},V,\beta } \right)}}{{{N_{\rm{e}}}!{N_{\rm{n}}}!}} $

      $ Q\left( {{N_{\rm{e}}},{N_{\rm{n}}},V,\beta } \right) = \sum\limits_\sigma {\int_V {{\rm{d}}q{\rm{d}}r\rho \left( {q,r,\sigma ,\beta } \right)} } $

      式中:q={q1, …,qNn}和r={r1, …,rNe}分别表示质子和电子的坐标;σ={σ1, …,σNe}表示电子的自旋态。

      密度矩阵可近似表示为

      $ \sum\limits_\sigma {\rho \left( {q,r,\sigma ,\beta } \right)} = \frac{1}{{\lambda _{\rm{n}}^{3{N_{\rm{n}}}}\lambda _\lambda ^{3{N_{\rm{e}}}}}}\sum\limits_{s = 0}^{{N_{\rm{e}}}} {{\rho _{\rm{s}}}\left( {q,\left[ r \right],\beta } \right)} $

      $ {\rho _{\rm{s}}}\left( {q,\left[ r \right],\beta } \right) = \frac{{C_{{N_{\rm{e}}}}^s}}{{{2^{{N_{\rm{e}}}}}}}\exp \left\{ { - \beta U\left( {q,\left[ r \right],\beta } \right)} \right\}\prod\limits_{l = 1}^n {\prod\limits_{k = 1}^{{N_{\rm{e}}}} {\varphi _{kk}^l\det {{\left| {\mathit{\Psi }_{ab}^{n,1}} \right|}_s}} } $

      电子的交换效应包含在(23)式最后的行列式中,因此密度矩阵有正有负,其绝对值用于MC的抽样计算,而符号则用于对物理量的计算中。Filinov等[116-118]利用该方法对氢(氘)的热力学性质进行了大量的计算。张其黎等[119]采用类似的方法研究了氢的物态方程。

      除了路径积分方法外,Pierleoni和Ceperley[120-121]发展了可用于有限温度模拟的CEIMC方法。利用BO近似可将(16)式表示为

      $ {Z_{{\rm{BO}}}} = \sum\limits_q {\oint {{\rm{D}}\mathit{\boldsymbol{R}}{\rm{exp}}\left\{ { - \int_0^\beta {{\rm{d}}t\left[ {{T_{\rm{n}}}\left( {\mathit{\boldsymbol{R}}\left( t \right)} \right) + {E_{\rm{q}}}\left( {\mathit{\boldsymbol{R}}\left( t \right)} \right)} \right]} } \right\}} } $

      $ {{\hat H}_{{\rm{el}}}}\left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{R}}} \right){\mathit{\Phi }_{\rm{q}}}\left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{R}}} \right) = {E_{\rm{q}}}\left( \mathit{\boldsymbol{R}} \right){\mathit{\Phi }_{\rm{q}}}\left( {\mathit{\boldsymbol{r}},\mathit{\boldsymbol{R}}} \right) $

      如果温度足够低,(24)式中可忽略激发电子能级,得到标准的BO近似形式,即离子在由电子基态能量E0(R)定义的势能面上运动。在高温下,必须对所有的电子能态求和,Cao等[122]将配分函数写成如下形式

      $ {Z_{{\rm{FEBO}}}} = \oint {{\rm{D}}\mathit{\boldsymbol{R}}{\rm{exp}}\left\{ { - \int_0^\beta {{\rm{d}}t\left[ {T\left( \mathit{\boldsymbol{R}} \right) + {E_{{\rm{el}}}}\left( {\mathit{\boldsymbol{R}}\left( t \right)} \right)} \right]} } \right\}} $

      (26)式被称为自由能BO近似(FEBO)。在BO近似下,将问题分解为两部分:电子部分完全包含在求势能面的问题中,可以由QMC方法求解;一旦求出势能面,则离子问题可由路径积分(量子离子)或经典(经典离子)MC方法求解。

      在CEIMC方法中,离子在有限温度下利用Metropolis MC算法求解,接受概率为

      $ A\left( {\mathit{\boldsymbol{R}} \to \mathit{\boldsymbol{R'}}} \right) = \min \left[ {1,\exp \left( { - \beta \Delta {E_{{\rm{BO}}}}} \right)} \right] $

      式中:R表示一组离子坐标;ΔEBO=EBO(R')-EBO(R),表示离子坐标改变后电子能量的变化。电子则采用基态QMC方法在零温下求解。最简单的基态QMC方法是VMC[123-124],VMC方法虽然不能得到足够精确的解,但可以为投影方法提供优化的试探波函数。在此基础上,可以采用DMC[125]或基态路径积分MC(Ground State Path-Integral, GSPIMC)[121]等投影方法进行更精确的计算。CEIMC方法已经用于对氢(氘)的研究[121, 126-128]

      不同的MC方法适用的热力学区域也不同。另外计算体系和模拟时间的不同也会导致模拟结果的差异,即所谓的收敛性。QMC方法用于氢氘物态方程的研究结果详见2.2节。

    • 为有效计算氢及同位素氘在较宽温度、压力区域的热力学性质并解决实际问题,需要较细致地考虑氢的离解、电离和绝缘体到金属的转变。解析模型可定性分析解决这些问题,但定量数据与精度较高的实验结果仍不一致。尽管当前数值模拟技术已有长足发展,氢的离解、电离和绝缘体到金属的转变仍是尚未解决的科学难题[129]。实用宽区域物态方程多采用半经验模型:按不同热力学区域的相结构进行建模并分别计算不同相的物态方程;基于化学图景,计算离化和电离系统的性质;利用实验数据进行定标;对化学组分不是很明确的区域,采用基于物理图景的QMD和QMC模拟获得理论数据;综合理论、实验和数值模拟结果,形成多模型集成的理论框架,理论框架的思想与金属材料的研究有类似之处[39, 47]

      相比较而言,涉及热力学区域较宽的代表性理论数据主要有3个。其一是SESAME数据,Kerley[130]给出的最早的数据是1972年版本,期间也发展了其他版本,据文献[39],2003年版本沿袭了早期的化学模型理论框架,改进了液体微扰论、分子振动和转动、离化平衡处理方法,采用现代密度泛函理论数值模拟结果和新的实验数据进行了定标,氘物态方程数据库的密度范围为0~100 g/cm3,温度范围为5~108 K,可用于ICF或其他工程设计,只对有限学术团体公开。其二是基于多模型集成的H-REOS.3数据库(以下简称REOS)物态方程数据,Becker等[48]给出了REOS物态方程数据表的理论模型示意图(见原文图 1),在离子强耦合和电子强简并的低温高密度和高温低密度区域,采用了QMD和QMC模拟数据,适用密度范围为3.0×10-8~1.8×103 g/cm3,温度范围为60~107 K,仅有氢的公开数据,主要目的是用于天体物理研究。从公开发表的资料看,REOS物态方程在部分区域给出的能量关系明显有问题。其三是SCvH数据,主要由Saumon等[131]、Ebeling[132]发展的化学模型计算大部分区域数据,同时结合TF模型,利用插值方法得到宽区域物态方程数据库,密度范围为10-6~3 000 g/cm3,温度范围为100~107 K,有公开数据,主要被用于天体物理(褐矮星、巨行星)模型研究。

    • 目前,QMD和QMC数值模拟是处理复杂多体系统最有效的方法,能得到很多实验无法测量区域(特别是离子强耦合和电子简并区域)的数据,是本文数据评估的基础。具体情况概述如下。

      Da Silva等[133]利用Nova激光加载实验获得氘的雨贡纽数据,其中:低压25 GPa的结果与气炮实验[28]及SESAME[39]、QEOS[134]等理论模型的结果一致;70~210 GPa高压数据明显偏软;150 GPa压力下的最大压缩密度ρmax>6ρ0,压缩比明显比理论预言结果大50%。为了从理论上认识液氘一次冲击压缩的极限压缩比是否能达到6,同时又受限于QMD模拟对计算机能力的要求,Lenosky等[135]首先给出了采用紧束缚近似的第一原理分子动力学模拟结果。随后Galli等[136]、Bagnier等[137]、Lenosky等[138]、Gygi等[139]、Desjarlais[140]、Bonev等[141]、Danel等[142]、Karasiev等[143-144]、Knudson等[12, 145]陆续给出了氘的物态方程QMD计算结果,Collins等[146]、Holst等[147]、Morales等[126]、Caillabet等[148]给出了氢的物态方程计算结果,刘海风等[149-150]计算了部分关心区域氘的物态方程,王聪等[151]发表了利用QMD和OFMD模拟的更宽区域的氢物态方程。这些研究表明,高温高压下,径向分布函数分析显示液氘系统确实存在分子离解,瞬时电荷密度等高线分析显示液体中存在原子、分子和多原子的分子团,电子态密度分析显示随着温度或密度的增加,液氘由分子绝缘体、半导体到金属流体转变。除Gygi等[139]采用CPMD方法的模拟结果外,其他QMD模拟得到的一次冲击压缩极限压缩比在4.5附近。上述数值模拟工作采用了不同的模拟程序,QMD方案也不尽相同,表 1详细列出了代表性工作的时间、第一作者、参考文献、方法、程序、机构和研究的温度、密度范围。

      Isotope Time Firstauthor Ref. Method Code Organization Temperature/K Density/(g·cm-3)
      D 1997 T. J. Lenosky [135] TBMD LANL 3 000-31 250 0.58-1.47
      D 2000 G. Galli [136] CPMD GP code LLNL 10 000 0.67,1.0
      D 2000 S. Bagnier [137] BOMD VASP CEA The effect of the local-spin-density-approximation functional
      D 2000 T. J. Lenosky [138] BOMD VASP LLNL,LANL 2 000-31 500 0.506-0.851
      D 2002 F. Gygi [139] CPMD GP code LLNL The compressibility is determined by shock-induced electronic excitations
      D 2003 M. P. Desjarlais [140] BOMD VASP SNL 3 800-50 105 0.553-1.756
      D 2004 S. A. Bonev [141] CPMD GP code LLNL 20-19 860 0.171-0.779
      D 2004 H. F.Liu [149] BOMD VASP IAPCM 1 000-10 000 0.506-1.0
      D 2016 J. F. Danel [142] QMD,OFWMD ABINIT,VAAQP CEA 11 604-116 045 0.2-20
      D 2016 V. V. Karasiev [143] QMD,OFMD PROFESS@Q-ESPRESSO,ABINIT Universityof Florida 2 000-250 000 0.2-10, No data list
      D 2017 M. D. Knudson [12] QMD VASP,Diff Exc SNL, WashingtonState University
      H 2001 L. A. Collins [146] BOMD VASP LANL,LLNL 5 000-30 000 0.334-0.525
      H 2008 B. Holst [147] BOMD VASP Institut für Physik,SNL 500-20 000 0.5-5.0
      H 2010 M. A. Morales [126] BOMD,QMC QBOX,CEIMC University of Illinois,University of L'Aquila 2 000-10 000 0.724-2.329
      H 2011 L. Caillabet [148] QMD,CEIMC PIMC, QMD, CEIMC CEA -116 045 0.2-5
      H 2013 C. Wang [151] QMD,OFMD ABINIT IAPCM 1.564×104-5.004441×107 9.82×10-3-1.347×103

      表 1  氢氘物态方程的QMD模拟工作

      Table 1.  Works of QMD simulation for equation of state of hydrogen and deuterium

      PIMC方法也已经广泛用于对氢(氘)的物性研究。1994年Pierleoni等[152]用RPIMC方法研究了高温高密度区(rs=0.5,1.0,1.61,2.0;T>5 000 K; 其中rs是以Bohr半径为单位的离子球半径)氢等离子体的物态方程,验证了简并参数大于0.4时Ichimarv解析物态方程的准确性,观测到低温低密度区分子的形成。1996年Magro等[112]用RPIMC方法研究了热稠密区(rs=1.75,1.86,2.0,2.2;T>5 000 K)氢分子的离解问题,定性支持化学模型的分子离解和原子离化过程。2000年至2001年,Militzer等[113-114]发展了变分密度矩阵节点近似,研究了氘密度分别为0.674和0.838 g/cm3、温度105 K≤T≤106 K的物态方程,评估了尺寸效应、积分时间步的影响,并计算了冲击雨贡纽,研究了低中密度(9.83×10-4 g/cm3ρ≤0.153 g/cm3,5 000 K≤T≤250 000 K)氢等离子体的热力学性质,分析了影响数值模拟的密度矩阵的准确性、分立时间误差、物态方程模拟结果的尺寸效应等,给出了与化学模型的定量比较。2011年胡素兴等[46]利用RPIMC方法计算了用于ICF模拟的氘的物态方程数据库(0.002 g/cm3ρ≤1 596 g/cm3,1.35 eV≤T≤5.5 keV)。

      2000年以来,Filinov等[116-118, 153-155]利用DPIMC方法对氢(氘)的热力学性质进行了计算,结合数值模拟结果,分析了特定温度和密度条件下离子和电子分布函数的变化特征,说明了可能存在氢等离子体相变,并对不同模型计算的热力学数据进行了交叉验证。张其黎等[119, 150]利用DPIMC方法研究了氢的物态方程。

      2004年, Pierleoni等[121]利用CEIMC方法研究了稠密金属氢(密度1.56~5.27 g/cm3,温度300~10 000 K)的热力学性质; 2010年, Morales等[126]计算了金属氢的物态方程(0.724~2.329 g/cm3,2 000~10 000 K);2015年, Tubman等[127]研究了沿Hugoniot曲线氘的分子-原子相变行为;李名锐等[128]利用CEIMC方法对液氘的一次冲击压缩特性进行了模拟。表 2详细列出了QMC数值代表性工作的时间、第一作者、参考文献、方法、程序、机构和研究的温度、密度范围。

      Isotope Time Firstauthor Ref. Method Organization Temperature/K Density/(g·cm-3)
      D 2000 B. Militzer [114] PIMC(VDM) University of Illinois at Urbana-Champaign 105-106 0.674,0.838
      D 2004 V. Bezkrovniy [155] DPIMC Universität Greifswald, Domstrasse, Germany;Russian Academy of Science 15 625-1 000 000 0.674,0.838,1.097
      D 2011 S. X. Hu [46] RPIMC University of Rochester,University of California 15 665-63 822 000 0.002-1 596
      H 20012003 V. S. Filinov [117]
      [116]
      DPIMC Russian Academy of Sciences 31 250-106 0.419
      H 2001 B. Militzer [114] RPIMC LLNL,University of Illinois at Urbana-Champaign 5 000-250 000 9.833×10-4-0.153
      H 2004 C. Pierleoni [121] CEIMC Universita' of L'Aquila, Via Vetoio, University of Illinois at Urbana-Champaign, Université Pierre et Marie Curie 300-10 000 5.267,2.697,1.561
      H 2010 M. A. Morales [126] CEIMC University of Illinois at Urbana-Champaign, University of L'Aquila, Italy 2 000-10 000 0.724-2.329
      H 2011 L. Caillabet [148] QMD,CEIMC CEA 116 045 0.2-5
      H 2015 Q. L. Zhang [150] DPIMC IAPCM 116 045 0.98×10-3-1 346.1

      表 2  氢氘物态方程的QMC数值模拟

      Table 2.  Works of QMC simulation for equation of state of hydrogen and deuterium

      根据1.2和1.3节的介绍,数值模拟方法本身和模拟预设参数等均会影响模拟结果。因此,本文首先对已公开发表的表 1表 2中大量的模拟数据进行对比分析,在2.1节和2.2节给出采用QMD和QMC数值模拟方法得到的物态方程的定量评估结果,结合对数值模拟结果的分析认识,在2.3节给出对几个典型宽区域模型物态方程与模拟数据相对差别的定量评估分析,最后得出对宽区域物态方程的定量可靠性的判断。

      表 1表 2的模拟数据大部分是一些离散的状态点,也有一些适用范围相对窄的解析拟合物态方程。如Holst2008[147]、Morales2010[126]、Caillabet的多相物态方程MP2011[148],其中:Holst2008物态方程直接由QMD数值模拟结果采用多项式拟合得到,适用密度范围0.5~5.0 g/cm3、温度范围500~20 000 K;Morales2010物态方程由CEMC的数值模拟结果拟合多项式得到,适用密度范围0.7~2.4 g/cm3、温度范围2 000~10 000 K;MP2011物态方程利用了第一原理分子动力学和耦合量子蒙特卡洛模拟结果,同时考虑了固相和液相结构的不同建模及物理极限,适用密度范围0.2~5.0 g/cm3、温度高至116 045 K。这些解析物态方程的共同特点是可以在拟合偏差范围内复现氢的数值模拟结果,可作为比较准确的模拟数据使用。在评估宽区物态方程时,常常缺少系统的数值模拟结果,利用这些解析物态方程可直接计算状态量,从而避免了由图形读取数据造成的偏差,或因数据点少而不能获得系统规律的问题。

      为了做定量比较,尽量选取包含明确的计算细节描述和状态量的文献数据。文献[148]指出氢氘物态方程可通过标度关系

      $ {\rho _{\rm{H}}} = \left( {1.008/2.014} \right){\rho _{\rm{D}}} $

      互推,但未见有严格证明。化学模型中关于氢氘振动和转动项之间并不能通过标度变换得到,同时数值模拟与核质量也有关系。因此本文对于具体数据点的应用尽可能忠实于文献,氢氘有所区别。仅在应用Holst2008、MP2011这两个物态方程评估氘的PIMC模拟数据时,采用该关系由氢的物态方程得到氘的物态方程。

    • Lenosky等[138]根据QMD模拟结果给出了氘的物态方程拟合函数结果,Holst等[147]给出了氢物态方程的拟合结果Holst2008,考虑到Lenosky的物态方程适用范围比较窄,选择Holst2008计算氢物态方程,利用(28)式的标度关系计算氘的物态方程,并与表 1中其他采用QMD方法模拟计算的氘物态方程结果进行对比,结果如图 1所示。由于各文献关于能量零点的选取不尽相同,这部分没有给出关于能量的比较结果。

      图  1  QMD模拟氘物态方程的压强相对偏差

      Figure 1.  Relative differences of pressure from QMD simulation for deuterium

      图 1可见,在Holst2008物态方程适用的密度0.5~5.0 g/cm3(对氘,1.0~10.0 g/cm3)、温度500~20 000 K范围内,Galli等[136]、Desjarlais[140]和Danel等[142]给出的物态方程的压强相对偏差较小,除个别点外,均小于10%。Danel采用QMD方法模拟了氘密度高达20 g/cm3、温度高达20 eV(1 eV=11 604.5 K)的物态方程,其高密度区的计算结果与Holst2008数据的相对偏差小于4%,说明Holst2008物态方程略作外推也还比较准确。相对偏差比较大的是Bonev和Lenosky的计算数据。Bonev等[141]计算结果的相对偏差最高可达35%。Bonev和Galli的计算采用了相同的GP(JEEP)程序,Galli严格考察了计算体系等的收敛性,其结果与Holst2008的相对偏差小于1%,由此我们判断Bonev为了做更大量的计算,选择的计算参数造成了比较大的偏差。图 1表明,温度2 000 K附近,Lenosky计算的压强与Holst2008物态方程的相对偏差均大于50%,其他几个温度的相对偏差均小于30%,大部分在20%以内。由于Holst和Lenosky均使用VASP程序,2 000 K附近标度关系的影响也不会太大,这个对比结果非常令人诧异,好在Knudson等[12]在2017年已明确指出Lenosky等[138]计算的压强并不收敛,可能是造成相对差异大的主要原因。图 1中Liu2004[149]是作者2004年的计算结果,对比可见压强的相对偏差约在10%以内,与作者当时选取了合适的收敛条件是一致的。

      同样地,以Holst2008的氢物态方程为标准,我们对比分析了Collins[146]、Morales等[126]和王聪等[151]计算给出的氢的物态方程,结果如图 2所示。可见虽然作者不同,采用的模拟程序也不尽相同,但计算得到的氢物态方程在温度1 000~30 000 K、密度0.3~5.0 g/cm3范围内的相对偏差仍然比较小,约在10%以内,这与模拟者本身对泛函、收敛条件等的把握是有关的。

      图  2  QMD模拟氢物态方程的压强相对偏差

      Figure 2.  Relative differences of pressure from QMD simulation for hydrogen

      以DFT为基础的QMD模拟方法中,电子的交换关联是通过近似方法处理的,关于氢的交换关联泛函的研究工作很多[69-79]。交换关联势的选取对物态方程计算结果的影响值得考虑。Knudson等[12]研究了不同密度泛函计算的氘的主雨贡纽曲线,发现在极限压缩度附近有微小的差别,在其附件材料中给出了基于不同泛函计算的氘的雨贡纽物态方程数据,这些数据与Holst2008物态方程的相对偏差如图 3所示。可见:取相同PBE泛函时,计算结果的相对偏差在5%以内;optB86b泛函的计算结果与PBE泛函类似;vdW-DF1和vdW-DF2两个泛函与PBE的计算结果在温度104 K以内、密度0.50~0.65 g/cm3,即偏离极限压缩比4.5(密度0.75 g/cm3)时,相对偏差最高可达28%,相对偏差随密度、温度的升高而降低。基于PBE泛函的QMD模拟得到的氘在20 GPa附近的雨贡纽压力比气炮实验结果偏低[138, 141],vdW泛函可能更好地考虑了氢分子的性质,与实验结果更接近[12],这与图 3的理论数据偏差情况是一致的,但目前缺乏基于vdW泛函的更宽热力学区域的数据。

      图  3  不同泛函计算的氘物态方程压强相对差别

      Figure 3.  Relative differences of pressure from QMD simulation for deuterium in different exchange-correlation analysis

      上述QMD方法中用到的交换关联泛函与温度无关,Karasiev等[143]采用QMD方法,选取有限温度局域密度泛函Karasiev-Sjostrom-Dufty-Trickey (KSDT)[144]和传统基态密度泛函LDA XC(Perdew-Zunger (PZ))两个不同的交换关联势,计算了氘物态方程,结果表明在研究的温度(低于2×105 K)、密度(0.2~10.0 g/cm3)范围内,不同交换关联势对压力的影响约为4%和6%,同时发现氢的主雨贡纽对显含温度的密度泛函不敏感。

    • 由于负符号问题,QMC方法多被用于研究高温物态方程,Holst2008物态方程的适用范围已不够用于这些物态方程的评估。选取MP2011为基准,适用于密度范围0.2~5.0 g/cm3(对氘,0.4~10.0 g/cm3)、温度高达116 045 K。图 4给出了胡素兴等[46]、Filinov等[117]、Militzer等[115]计算的氘物态方程数据的相对偏差分析。图 4表明,在10 000 K时,Militzer采用RPIMC的计算结果相对偏差可达50%左右,随着温度升高,其相对偏差缩小到10%以内。在15 625 K时,Filinov采用DPIMC的计算结果相对偏差可达92%,温度高于60 000 K时,该相对偏差下降到10%以内。温度低于60 000 K时,胡素兴采用RPIMC的计算结果相对偏差可高达30%。张其黎基于DPIMC的研究表明,在大部分区域压强与胡素兴RPIMC结果的偏差在10%以内,温度较低时有部分偏差超过20%;同时指出了RPIMC两个压强数据点的错误。在离解和电离区,DPIMC压强位于TF和TFC模型之间,在离解区与RPIMC的相对偏差达到15%,在电离区两者一致。在低密度区,DPIMC的内能与RPIMC的内能接近,小于QMD和IG模型;在高密度区,DPIMC的压强位于TFC和IG模型之间,内能则略小于IG模型。PIMC模拟结果在低温下表现出明显的分散性,原因有两个。一是为了解决负符号问题,Militzer、胡素兴采用了变分密度矩阵节点近似,该方法利用固定节点近似避免了负符号问题,但引入了两个不可控因素:(1)固定节点限制了对电子交换的准确描述, (2)节点本身是近似的。二是在量子效应比较明显的区域,这些因素对模拟结果会有较大的影响。胡素兴模拟的最低密度为0.506 g/cm3,最高密度为10.0 g/cm3,在温度高于MP2011适用范围的116 045 K区域,其模拟结果与MP2011的相对偏差约在10%以内,原因是Caillabet在做MP2011物态方程时,考虑了高压高温极限条件,因此具有较好的外推性能。

      图  4  PIMC模拟氘物态方程的压强相对偏差

      Figure 4.  Relative differences of pressure from PIMC simulation for deuterium

      Morales等[126]和Pierleoni等[121]利用CEIMC、Filinov等[117]利用DPIMC方法计算了氢的物态方程,图 5给出了其结果与MP2011物态方程的相对偏差。由图 5可见,DPIMC方法计算的氢的物态方程在温度31 250 K时,相对偏差高达44%,与前面氘的情况类似。Morales计算结果的个别点的相对偏差大于10%,其他在10%以内。Pierleoni在两个低密度1.561和2.697 g/cm3下的计算结果相对偏差小于10%,在高密度5.267 g/cm3、温度低于10 000 K时,6个点的压强相对偏差均大于35%。Morales和Pierleoni的模拟中关于零温电子基态的处理方法不同,分别采用了RMC、VMC方法,二者计算的物态方程在10 000 K以下、密度较低时结果基本一致。图 1显示Holst2008在高密度5~10 g/cm3、数万开尔文附近的外推结果与QMD模拟结果一致,在此密度、温度范围内,MP2011与Holst2008物态方程接近,也就是说MP2011物态方程具有较高的置信度,图 5中高密度5.267 g/cm3下Pierleoni等[121]基于VMC计算的物态方程结果与MP2011物态方程的偏离是难以理解的。

      图  5  PIMC模拟氢物态方程的压强相对偏差

      Figure 5.  Relative differences of pressure from PIMC simulation for hydrogen

      第一性原理计算常被认为是求解过程没有可调经验参数的方法。2.1节与2.2节的比较及1.1节和1.2节的方法介绍表明,由于实现第一性原理计算的理论方案仍有不可避免的近似,以及程序使用者对物理方法的理解和初始设置参数的选择不同,需要仔细判断计算结果的可靠性。

    • REOS[48]是目前可以获得的适用温度、密度区域最宽的氢的物态方程。本节结合MP2011物态方程[148]、QMC[46]和DFT(QMD+OFMD)[151]模拟结果及改进的自由能模型对其进行定量分析。

      前两节已对多个研究小组采用QMD和QMC数值模拟方法得到的氢氘物态方程压强数据进行了系统、定量的对比分析,结果表明MP2011物态方程与大多数采用QMD和QMC的数值模拟结果一致。进一步比较了Holst2008、Morales2010和MP2011物态方程的能量,发现MP2011物态方程计算的能量与Holst2008的能量偏差在分子离解区域高于15%,某些状态点的相对差异高达几百倍以上,但与Morales2010给出的能量的相对偏差小于5%,表明MP2011物态方程能较好地拟合CEIMC模拟得到的分子离解区的能量关系。

      在MP2011物态方程适用的温度密度范围内,图 6给出了MP2011与REOS的定量比较。图 6表明:压强相对偏差在大部分区域小于10%,在两个蓝色区域为10%~15%之间;在104 K以内、密度1~3 g/cm3的能量相对偏差也在10%以内,没有出现类似Holst2008物态方程所呈现的明显的分子离解导致的能量变化特征;低密度、105 K温度附近的能量相对偏差变化剧烈,变化区域随着密度的增大温度逐渐降低,结合PIMC方法模拟得到的分布函数及自由能模型分析,此区域的能量变化关系源自原子电离。

      图  6  氢的MP2011与REOS物态方程的相对偏差

      Figure 6.  Relative differences of equation of state for hydrogen between MP2011 and REOS

      胡素兴等[46]采用RPIMC方法系统计算了氘的物态方程,本文采用标度关系由其计算氢的数据,并与REOS物态方程进行定量比较,结果见图 7。由于RPIMC和REOS采用了不同的能量零点,能量的比较结果在相对差异上加0.5得到。由图 7可见:在T>105 K的大部分区域,二者基本一致,压强和能量相对偏差在2%以内;T≤105 K, REOS与RPIMC模拟结果的能量相对差异在个别区域可高达20%以上;在温度104 K、密度1 g/cm3附近,RPIMC模拟压强及能量的相对偏差在10%~20%之间。高温高密度的深蓝色区域的差异是因REOS数据范围(最高密度、温度分别为1.8×103 g/cm3和1.0×107 K)较窄,外插引起的偏差所致。

      图  7  RPIMC模拟的氘物态方程与REOS物态方程的相对偏差

      Figure 7.  Relative differences of equation of state for deuterium between RPIMC and REOS

      王聪等[151]的DFT模拟结果和REOS采用了不同的能量零点,将DFT能量减去15.502 eV/atom进行对比。图 8给出了DFT模拟结果与REOS的定量比较,可见:在温度T>105 K的大部分区域,压强相对偏差在2%以内;在温度105 K附近及以下,压强的相对偏差在2%~8%之间;在T>105 K的区域,能量的相对偏差在低密度处大多小于2%;在高密度及高温的大部分区域,能量相对偏差在2%~8%之间。这和上面RPIMC的比较结果不同,分析是由于DFT模拟采用了QMD和OFMD两种模拟方法,二者能量零点并不相同所致。在温度105 K附近及以下,部分区域能量的相对偏差高达20%以上。同样地,高温高密度的深蓝色区域的差异是因REOS数据外插所致。

      图  8  DFT模拟的氢物态方程与REOS物态方程的相对偏差

      Figure 8.  Relative differences of equation of state for hydrogen between DFT and REOS

      基于改进的自由能模型,结合插值方法,也得到了氢的宽区域物态方程。计算结果与REOS的对比如图 9所示。Becker等[48]给出了REOS物态方程数据表的理论模型示意(见原文图 1),结合图 9相对偏差的分布情况可见,改进的自由能模型可统一代替FVT、SCvH、CP模型,在这些模型适用的区域,与REOS压强及能量的相对偏差最大约为8%,大部分区域的差别在2%以内。REOS采用DFT-MD构建密度0.1~100 g/cm3、温度低于105 K区域的物态方程,改进的自由能模型采用Pade近似描述等离子体的库仑相互作用,计算给出的10~100 g/cm3区域的计算结果与REOS结果基本一致,即与QMD一致,相对偏差小于2%。

      图  9  改进的自由能模型计算的氢物态方程(LiEOS)与REOS物态方程的相对偏差

      Figure 9.  Relative differences of equation of state for deuterium between modified free energy model (LiEOS) and REOS

    • 针对QMD、QMC数值模拟方法和自由能模型用于氢氘物态方程研究中的普遍性问题,文献[11]已有较多点评,这里将结合2.2节~2.4节的定量评估结果,进一步做具体讨论。

      QMD和QMC数值方法将电子和核当作基本粒子,它们之间通过库仑势相互作用,用量子理论描述原子和分子的形成及其统计性质,通过大规模的数值模拟确定材料的性质,理论上而言,这两种方法是严格精确的。

      实际应用层面,QMD方法用于氢氘的物态方程模拟,也获得了一些很好的结果[136, 140, 142, 147, 149, 151],但也存在一些固有的缺陷。首先,QMD模拟中电子的描述基于DFT,其中交换关联泛函是最不确定因素[11-12, 71-79],需要通过实验和其他更严格的方法验证其准确性,对于氢氘而言,为获得相空间中不同热力学状态的准确描述,需要多个交换关联函数,如对低压(冲击加载低于20 GPa)区域,vdW泛函比PBE泛函计算的压强高约20%以上,如图 3低密度区所示,与实验结果的对比表明vdW泛函比PBE描述更准确。其次,常用的基于轨道方法的BOMD方法在电子结构的自洽计算中,需要求解与单粒子轨道相关的KS方程,随着温度的升高,需要计算的轨道数目迅速增加,导致计算量大大增加,因此,轨道方法并不适合很高温度下的氢氘物态方程模拟。王聪等[151]为解决这个问题,在高温高密度区采用无轨道的OFMD方法作为补充,获得了图 8(b)“+”号所示温度、密度状态点的物态方程,但又引入了新的问题,即OFMD和QMD能量零点不统一,由此带来的能量系统差别在2%~8%之间。另外,因为交换关联泛函的近似、氢高压结构的复杂和不确定性,QMD方法模拟氢的金属化、等离子体相变(PPT)等问题仍是开放的课题[129, 156],这些问题的进一步深入研究对构建氢的宽区域多相物态方程是十分重要的。

      QMC方法用于模拟电子-离子系统时,电子和离子都用路径积分描述,可自动处理质子的零点运动,理论上可以不采用BO近似,由于不存在单粒子近似,因此对电子交换关联的描述是准确的。与QMD方法相比,其最大的优点是可以准确处理分子的离解,MP2011[148]、Holst2008[147]和CEIMC[126]3个物态方程的能量对比表明,分子离解的准确描述引起的能量变化在10%以上。RPIMC[46]用于氘的物态方程模拟,在高温高密度区获得了理想的结果,如图 9(b)“+”号所示温度T>105 K区域的物态方程。但在T<105 K区域,RPIMC[46, 115]、DPIMC[117]及CEIMC[121, 126]模拟结果仍存在较大差异,如图 4图 5所示,影响主要来源于费米符号问题,即交换的反对称性导致积分的正负项相互抵消。为解决这个问题,上述3个蒙特卡洛方法分别引入新的近似,但这些近似方法对计算精度的影响有多大,仍是一个开放的课题。另外,由于质子和电子质量相差3个量级,在相空间的运动路径相差甚远,因此电子、离子的全蒙特卡洛模拟存在统计抽样困难,准确的计算对计算资源的要求非常高。

      改进自由能模型很好地给出了氢物态方程的低密度、高温和高密度极限计算结果,如图 9所示。另一方面,自由能模型将分子和原子作为系统的基本成分,实际上是把电子和离子的束缚态作了近似处理,在特定的温度、密度下,原子和分子的寿命随温度、密度的增加衰减得很快,且电子和离子的平均距离变得与键长可以比拟,因此原子和分子的定义将变得不清楚,这类方法的准确性需要通过其他独立的数值模拟方法进行检验。在低密度(9.833×10-4 g/cm3ρ<0.153 g/cm3)、高温(5 000 K<T<250 000 K)[114]和高温(T>105 K)[21]区域,RPIMC[46]模拟结果检验了基于自由能模型计算的氢物态方程。QMD计算检验了温度T<105 K、密度10~100 g/cm3[48]区域自由能模型的计算结果。

      REOS的氢物态方程与RPIMC、QMD和改进自由能模型的结果对比见图 7~图 9。分析认为REOS数据存在如下问题:高温高密度区域数据范围不够,可导致外插偏差;温度103~105 K之间能量变化有明显差异,涉及分子离解和原子离化,需要采用统一理论框架重新处理;温度103~104 K之间、密度1.0 g/cm3左右需采用更准确的QMD数据重新定标;低密度区能量计算不合理。另外,REOS只有氢的物态方程,从改进的自由能模型看,分子的振动和转动对物态方程的贡献不满足(28)式的标度变换关系。氢和氘的宽区物态方程需要分别构建。

      上述对比分析结果为采用多分区和多模型优化集成方式构建氢氘宽区域物态方程奠定了理论和数据基础,由于篇幅关系,关于氢氘宽区域物态方程的构建和验证将结合改进的自由能模型另文讨论。

    • 根据对多个研究组不同时期的QMD、QMC和自由能模型计算的氢氘物态方程进行定量分析评估,得到如下结论。

      (1) 第一原理是计算物态方程的重要工具,考虑到模拟方法的理论基础,密度在零点几到几克每立方厘米时,QMD适用于温度从几千到几十万开尔文的区域,QMC适用于温度高于几万开尔文,CEIMC适用于几百到万开尔文,具体适用温度随密度而变化。结合多种数值模拟方法,如果能承受昂贵的计算机,理论上可得到氢氘的宽区域物态方程。

      (2) 实际操作方面,受限于目前的计算机硬件及模拟能力,仍然只能采用各种不同层次的近似,对有限的热力学相空间进行模拟研究。在模拟初始参数控制较好的情况下,得到的结果有些与实验一致,有些与实验仍有较大的差别。特别是在实验不确定度较小时,其可为理论方法提供更好的检验。

      (3) 目前,基于第一原理数值模拟得到的氢、氘物态方程只能描述有限的热力学相空间,同时,由于不可避免的物理或数值近似,还无法替代半经验模型的宽区域物态方程。基于多模型集成的REOS数据库,在105 K以下与模拟结果的相对偏差仍比较大,数据点不够密,最高温度仅到107 K,由REOS数据库不能直接得到氘的数据库。氢、氘物态方程均不能满足工程应用。需采用多模型优化集成的技术路线,充分发挥解析模型和数值模拟的优点,结合高精度的实验研究,分别构建数据库形式的氢、氘宽区物态方程。

参考文献 (156)

目录

    /

    返回文章
    返回