为什么我尝试计算的指数即使在我的 Fortran 程序中使用双精度也会趋于零？

Question

我正在尝试编写一个程序，使用热力学配分函数计算反应的动力学常数。然而，当尝试计算指数时，即使使用双精度和四精度，指数也太小并且会趋于零。现在我不知道这是概念错误还是编程错误，但无论如何我都想要第二个意见。

我附上了我正在使用的程序和输入文件，如果格式很尴尬，就像我第一次发帖一样，很抱歉。

输入文件

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.43e10，但我无法继续计算指数。 欢迎任何建议，如果需要，愿意提供更多信息。

Answer 1

通常问题出在单位和转换上。你应该检查一下他们是否正确。另一个问题可能是您编写指数的方式，您可以尝试使用“E”格式。 real*16 应该足以执行此类计算。如果这还不够，您应该检查是否有任何库可以提高精度。 我希望这能有所帮助，让我们保持最新状态