双/四精度计算指数仍归零？Eyring反应速率常数求解求助

算力阿峰

2026-6-24

Eyring方程分子版本速率常数计算指数项下溢问题

我正在编写程序，通过热力学配分函数结合Eyring方程计算反应速率常数。摩尔版本计算正常，但分子版本中指数项数值过小导致结果直接归零，即使使用double（real*8）和quadruple（real*16）精度也无法解决。这不是精度不足的问题，而是指数函数的特性导致，需要能获得有效数值结果的方案。预期分子版本结果约为1.43e10。

输入文件

298 !temperatura, K
6570 !delta H formazione, J/mol

1 
18.9984  !massa rg1 au
4 ! deg1

2
2.016 !massa rg2 au
4395.2 !vibr, cm-1
0.277 !inerzia au**2
1 !deg

3
linear
21.014 !massa tr au
4007.6 397.9 397.9 310.8  !vibr cm-1
7.433
4 !deg
2 !numero simm

程序代码

program cinetica

c DICHIARAZIONE COSTANTI

      real*8 ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 T, deltaf, qpart(3), q, k2
      real*8 part1, part2
      real*16 part3, esp
      common /cost/ ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      integer atom

      ktransm= 1 !coefficiente trasmissione
      kboltz= 1.38*10**(-23.) !costante di boltzmann, J/K
      rgas= 8.31 ! costante universale dei gas, J/mol*K
      hplanck= 6.62*10**(-34.) !costante di planck, J*s
      AvN= 6.022*10**(23.) ! costante di avogadro, 1/mol
      Pst= 10**(5.) ! pressione standard, Pa
      c= 3.*10**(10.) ! velocità della luce in cm/s
      
C APERTURA FILE E LETTURA DATI
      open(5, file="dati.in")
      open(10, file="cinetica.out")

      read(5,*) T  !temperatura in K
      read(5,*) deltaf  !delta H di formazione in J/mol

      Vst= (rgas*T)/Pst
      pi= 2*asin(1.)

c legge masse con un do e in base al numero di atomi sceglie come calcolare numero di partizione, alla fine produco q di ognuno che assegna a un vettore q con tutte le funzioni
      do 50 i=1,3
        read(5,*) atom    !legge numero di atomi
        write(6,*) atom

        if(atom.eq.1) then
          call monoatomic(T,q)

        elseif(atom.eq.2) then
          call diatomic(T,q)

        else
          call polyatomic(atom,T,q)

        endif

        qpart(i)=q

 50   continue

      write(6,*) "funzioni part", (qpart(n), n=1,3)


C CALCOLO COSTANTE CINETICA (EQUAZIONE DI EYRING) CON FUNZ. PARTIZIONE DA SUBROUTINE E COSTANTI DICHIARATE IN PRECENDENZA
  
c   versione molare, funziona?

c      part1= (ktransm*kboltz*T*Vst)/hplanck
c      part2= (AvN*qpart(3))/(qpart(1)*qpart(2))
c      part3= exp(-deltaf/(rgas*T))
  
c      k2= part1*part2*part3
c      write(6,*) "costante cinetica molare", k2

c   versione molecolare, esponente va a zero

      part1= (kboltz*T)/hplanck
      part2= (qpart(3))/(qpart(1)*qpart(2))
      esp=-deltaf/(kboltz*T)
      part3= exp(esp)

      write(6,*) "esponente e esponenziale", esp, part3

      k2= part1*part2*part3
      write(6,*) "costante cinetica molecolare", k2
      
      stop
      end
c-----------------------------------------------------
C SUBROUTINE MONOATOMICO

      subroutine monoatomic(T,qtot)
      common /cost/ ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 T, m, deg, qtr, qel, qtot

      read(5,*) m
      read(5,*) deg
      
      write(6,*) m, deg
      m=m*(1.67*10**(-27.))
      qtr= ((2*pi*m*kboltz*T)/(hplanck**(2.)))**(3./2.)*Vst
      qel=deg

      qtot=qtr*qel
      write(6,*) qtr, qel, qtot


      return
      end



c-------------------------------------------------------------
C SUBROUTINE DIATOMICO
      subroutine diatomic(T,qtot)
      common /cost/ ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 T, m, v, deg, qtr, qrot, qvib, qel, qtot
      real*8 i

      read(5,*) m
      read(5,*) v
      read(5,*) i
      read(5,*) deg

c conversione unità di misura per massa e inerzia
      m=m*(1.67*10**(-27.))
      i=i*(1.67*10**(-27.))*(1.67*10**(-27.))


      qtr= ((2*pi*m*kboltz*T)/(hplanck**(2.)))**(1.5)*Vst
c      write(6,*) "Vst", Vst
      qrot=(8*pi**(2.)*i*kboltz*T)/(hplanck**(2.))
      qvib=1/(1-exp((-hplanck*c*v)/(kboltz*T)))
      qel=deg
      
      qtot=qtr*qrot*qvib*qel
      write(6,*) qtr, qrot, qvib, qel, qtot

      return
      end


c----------------------------------------------------------------------
C SUBROUTINE POLIATOMICO

      subroutine polyatomic(atom,T,qtot)
      common /cost/ ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 ktransm, kboltz, rgas, hplanck, AvN, Pst, pi, Vst, c
      real*8 T, m, v(500), deg, qtr, qtot, qvib, qvib1, qel
      real*8 i
      character*20 geom
      integer modes, symm, atom

      read(5,*) geom

      if(geom.eq."linear") then
         modes=3*atom-5
      elseif(geom.eq."non linear") then
         modes=3*atom-6
      else
         write(6,*) "check input: define geometry (linear/non linear)"
      endif

      read(5,*) m
      read(5,*) (v(n), n=1, modes)
      read(5,*) i
      read(5,*) deg
      read(5,*) symm

c conversione unità di misura per massa e inerzia
      m=m*(1.67*10**(-27.))
      i=i*(1.67*10**(-27.))*(1.67*10**(-27.))

      qtr= ((2*pi*m*kboltz*T)/(hplanck**(2.)))**(3./2.)*Vst
      qrot=(8*pi**(2.)*i*kboltz*T)/(symm*(hplanck**(2.)))
      qel=deg

      qvib=1

      do 60 nn=1, modes

        qvib1=1/(1-exp((-hplanck*c*v(nn))/(kboltz*T)))   !parziale del calcolo vibrazionale
        qvib=qvib*qvib1     !produttoria

 60   continue

      qtot=qtr*qrot*qvib*qel
      write(6,*) qtr, qrot, qvib, qel, qtot

      return
      end

输出结果

1
   18.998400000000000        4.0000000000000000     
   2.0035785894328219E+030   4.0000000000000000        8.0143143577312874E+030
           2
   6.9257471645836213E+028   5.7237652847843724E-007   1.0000000006050396        1.0000000000000000        3.9641351215822116E+022
           3
   2.3307405557255242E+030   7.67955771E-06   1.7660399314814128        4.0000000000000000        1.2644179482433861E+026
 funzioni part   8.0143143577312874E+030   3.9641351215822116E+022   1.2644179482433861E+026
 esponente e esponenziale  -1597607262508866292678656.00000000000         0.00000000000000000000000000000000000      
 costante cinetica molecolare   0.0000000000000000

问题分析与解决方案

核心错误

当前代码中计算分子版本指数项时，直接使用了摩尔焓变（deltaf=6570 J/mol）除以单个分子的玻尔兹曼常数与温度的乘积，导致指数项数值为$-1.6e24$，远超出浮点数的可表示范围（double最小指数约为-709，quadruple约为-16383），最终exp(esp)下溢为0。

正确的做法是将摩尔焓变转换为单个分子的能量：$\Delta E = \Delta H / N_A$，再代入指数项计算。

具体解决方案

1. 修正指数项计算（最直接）

修改分子版本的指数项代码：

esp= - (deltaf / AvN) / (kboltz*T)

修正后指数项数值约为$-2.65$，处于浮点数有效范围内，exp(esp)不会下溢，计算结果可得到预期的1.43e10。

2. 利用摩尔版本结果推导

已知摩尔版本与分子版本速率常数的关系：$k_{molec} = k_{mol} / N_A$。先计算正常运行的摩尔版本结果，再除以阿伏伽德罗常数得到分子版本结果，完全规避指数项下溢问题。

3. 对数计算法（通用数值稳定方案）

将所有项转换为对数形式计算，再取指数还原：

log_k = log((kboltz*T)/hplanck) + log(qpart(3)) - log(qpart(1)) - log(qpart(2)) + (-deltaf/AvN)/(kboltz*T)
k2 = exp(log_k)

这种方式避免了大数值与极小数值直接相乘，提升计算稳定性。

内容的提问来源于stack exchange，提问作者Giulia Gamboni

双/四精度计算指数仍归零？Eyring反应速率常数求解求助

输入文件 #

程序代码 #

输出结果 #

问题分析与解决方案 #

核心错误 #

具体解决方案 #