TD calculations using ARIR

Ther­mo­dy­nam­ic cal­cu­la­tions tak­ing into account the inter­nal rota­tions.
1. An auto­mat­ed pro­ce­dure

MOLTRAN can per­form the ther­mo­dy­nam­ic cal­cu­la­tions tak­ing into account the inter­nal rota­tions con­tri­bu­tions. It can be done or in a man­u­al mode (using the com­mands in  .tdi-file) either in com­plete­ly auto­mat­ic mode. The auto­mat­ic mode is much more sim­ple and use­ful. How­ev­er, it will prob­a­bly require some addi­tion­al inspec­tion of the results in order to ver­i­fy some pro­gram deci­sions. Below, the results of the auto­mat­ic mode cal­cu­la­tions are reviewed using the SiF3-OH mol­e­cule as an exam­ple. 

To obtain the pre­sent­ed results, use the above-men­tioned sam­ple file Moltran-Sample1.txt. Start MOLTRAN and open the file Moltran-Sample1.txt using the open­ing dia­log or run the com­mand from a com­mand string: 

moltran Moltran-Sample1.txt

Choose menu items: “Calculations”->“Thermodynamics”. In the win­dow appeared check “Exit pro­gram after TD cal­cu­la­tions” and click “OK”. The pro­gram will fin­ish the work.

You can also use the way not requir­ing the inter­ac­tive oper­a­tions. Just run the com­mand:

moltran Moltran-Sample1.txt /0 /arir

After the pro­gram is fin­ished, the out­put file (Moltran-Samle1.log) con­tains the results of ther­mo­dy­nam­ic cal­cu­la­tions tak­ing into account the con­tri­bu­tion of sin­gle inter­nal rota­tion of OH group in SiF3OH. 

The cal­cu­la­tion is per­formed auto­mat­i­cal­ly, using  the Auto­mat­ic Recog­ni­tion of Inter­nal Rota­tions (ARIR) mod­ule imple­ment­ed in MOLTRAN. The ARIR pro­ce­dure uses an orig­i­nal algo­rithm. The method for TD account­ing of inter­nal rota­tions based on the the­o­ry devel­oped by Pitzer, Gwinn, and Kil­patrick. (see J.Chem.Phys.1942,10).

Describ­ing the out­put, we skip the begin­ning of file and the stan­dard ther­mo­chem­istry data and begin the dis­cus­sion from the results of inter­nal rota­tion treat­ment. The data of inter­est are start­ed with a head­er:

*** Auto­mat­ic Recog­ni­tion of Inter­nal Rota­tions ***

The ARIR pro­ce­dure starts at this point. At first, it tries to iden­ti­fy “basis rota­tions”, i.e. rota­tions which do not change the con­nec­tiv­i­ty matrix. In the cur­rent exam­ple, there is the only can­di­date cor­re­spond­ing to the rota­tion around atoms 1 and 5:

 *** Automatic Recognition of Internal Rotations ***

 Possible Internal Rotations :  1
 Int.Rot.  1   Axis:  1  5   NRot:  1   Atoms:  2  Group:  5  6

 Atom  XYZ   Internal Rotations
                1
 Si1    X    0.0000
        Y   -0.0001
        Z   -0.0112
 F2     X    0.0001
        Y    0.0000
        Z    0.0084
 F3     X    0.0237
        Y   -0.0074
        Z   -0.0053
 F4     X   -0.0238
        Y    0.0074
        Z   -0.0055
 O5     X    0.0001
        Y   -0.0003
        Z   -0.0405
 H6     X   -0.0021
        Y    0.0088
        Z    0.9984

Then, the nor­mal modes are expand­ed in the basis rota­tions. Two of them (vibra­tional modes 7 and 9 with fre­quen­cies 113 and 276 cm-1) has expan­sion coef­fi­cients exceed­ing the pre­de­fined thresh­old:

  Vibrational modes expanded to internal rotations. (Int.Rot.threshold = 0.50)
  Mode  Frequency       Internal Rotation Coefficient
   No.     cm-1            1
    1      0.00         0.0000            
    2      0.00         0.0000            
    3      0.00         0.0000            
    4      0.00         0.0000            
    5      0.00         0.0000            
    6      0.00         0.0000            
    7    113.42        -0.9961

From these two modes MOLTRAN tries to choose only one mode (because there is only one basis rota­tion) which is most sim­i­lar to a “pure inter­nal rota­tion”. This pro­ce­dure is inter­nal­ly ambigu­ous, even if we try to do that man­u­al­ly, because the “pure rota­tion” is a kind of ide­al­iza­tion and we can not to rep­re­sent any mol­e­c­u­lar mov­ing by rota­tion only. The results of such assign­ment should be accept­ed with a cau­tion. It is rec­om­mend­ed to test the assign­ment using the MOLTRAN vibra­tion ani­ma­tion pro­ce­dure. Name­ly, check, whether the vibra­tion is real­ly sim­i­lar to the vibra­tion of one part of mol­e­cule rel­a­tive­ly to anoth­er one.

 Vibrational modes which are most similar to internal rotations:
 (please check the program decision manually)
   No. Mode    Freq.         Internal Rotation Contributions (normalized)
                               1
    1    7    113.42        -1.0000

After the rota­tion mode is deter­mined, MOLTRAN esti­mates the rota­tion bar­ri­er height. Using ARIR, the semi­em­pir­i­cal pro­ce­dure is uti­lized (PM3 by default), and the results are print­ed below:

 Estimation of internal rotation potential by PM3 method

 Relative potential energies in kJ/mol as a function of rotation angle
 No.  Angle   Internal rotation
      (deg)       1
  1    0.00    0.000000
  2   10.00    0.045191
  3   20.00    0.178237
  4   30.00    0.383180
  5   40.00    0.634619
  6   50.00    0.901563
  7   60.00    1.154014
  8   70.00    1.370036
  9   80.00    1.540072
 10   90.00    1.666867
 11  100.00    1.762051
 12  110.00    1.841670
 13  120.00    1.921974
 14  130.00    2.015356
 15  140.00    2.126254
 16  150.00    2.247998
 17  160.00    2.362756
 18  170.00    2.446492
 19  180.00    2.478258
 20  190.00    2.449524
 21  200.00    2.368188
 22  210.00    2.254823
 23  220.00    2.133464
 24  230.00    2.022228
 25  240.00    1.928187
 26  250.00    1.847272
 27  260.00    1.767353
 28  270.00    1.672318
 29  280.00    1.546128
 30  290.00    1.377003
 31  300.00    1.161891
 32  310.00    0.909959
 33  320.00    0.642805
 34  330.00    0.390276
 35  340.00    0.183442
 36  350.00    0.047936
 37  360.00    0.000000

Now, the inter­nal rota­tion con­tri­bu­tion is esti­mat­ed. The mol­e­c­u­lar para­me­ters used for this esti­ma­tion is print­ed as fol­lows:

 *** Contributions of internal rotations ***

 *** Internal rotation No.  1 (replacing the vibration mode  7) 

 Masses (amu) and moments of inertia (amu*bohr^2) of rotating groups:
 Group 1:  M1=  17.00274   I1=    2.45338
 Group 2:  M2=  84.97213   I2=  442.72681
 Reduced I=I1*I2/(I1+I2)     =    2.43986
 Pitzer's estimator I0       =    2.32452
 Pitzer's estimator I        =    2.32452
 Rotation axes:
 Point 1:  Atom  5  Coords:  3.039279  0.248852  0.004620
 Point 2:  Atom  1  Coords:  0.005754  0.014445  0.000036

The cal­cu­lat­ed rota­tion poten­tial is expand­ed in a Fouri­er series in order to deter­mine its fold­er­ness (rota­tion sym­me­try) to der­mine the sym­me­try num­ber need­ed in ther­mo­dy­nam­ic cal­cu­la­tion. This deter­mi­na­tion is based on the cal­cu­la­tion of max­i­mum val­ue of expan­sion coef­fi­cients in the Fouri­er series. Unfor­tu­nate­ly, it is also very ambigu­ous pro­ce­dure because the true poten­tial is fre­quent­ly masked by addi­tion­al inter­ac­tions. In this exam­ple, the pro­ce­dure is failed to get cor­rect sym­me­try of rota­tion poten­tial: the obtained val­ue is 1 instead of 3 (3-fold poten­tial of rota­tion). This is a result of sig­nif­i­cant asym­me­try of poten­tial obtained at the PM3 lev­el. It is inter­est­ing that the ab ini­tio poten­tials are fre­quent­ly more sym­met­ric and the pro­ce­dure gives the cor­rect val­ue of sym­me­try num­ber. The sym­me­try num­ber can be also be set man­u­al­ly using the TDI-file as described below.

 Rotation barrier is estimated from the Fourier-transformed rotational potential
 V = c0/2 + c1*Cos(x) + c2*Cos(2*x) + ... ,    0<=x<=Pi
 First coefficients are:
     n      C(n)
     0    25.83821
     1    -9.71418

The rota­tion bar­ri­er height is esti­mat­ed on the basis of max­i­mum val­ue of poten­tial curve. Because the obtained bar­ri­er height is low­er than 1.4RT at the giv­en tem­per­a­ture, the con­clu­sion is made that this rota­tion is free:

 Rotation barrier height for n-fold potential V=V0/2*[1-Cos(n*x)] (n=1):
 V0 =     2.4783   kJ/mol
          0.5923 kcal/mol
          0.9997 RT   (T=  298.15K) => Free rotation (V0<1.4RT)

Now, the ther­mo­dy­nam­ic con­tri­bu­tion is esti­mat­ed. Free-rota­tion par­ti­tion func­tion used in Pitzer-Gwinn tables is print­ed out. Then the con­tri­bu­tion of the cur­rent rota­tion mode is giv­en in dif­fer­ent approx­i­ma­tions. The most com­pre­hen­sive is the last col­umn (obtained on the basis of Pitzer-Gwinn’s numer­ic tables for com­pre­hen­sive esti­ma­tor of iner­tia moments).

 Free-rotation function: 1/Qf =     1.946788E-01   Qf =     5.136666E+00
 Free-rotation (Pitzer): 1/Qf =     1.994503E-01   Qf =     5.013780E+00


 Thermodynamic functions for internal rotation  1 (vibration mode  7)

                  Harmonic   Free Rotator   Pitzer-Gwinn tables 
                 Ocsillator  I(*)    I(**)     I(*)    I(**)
 C,    J/K/mol=    8.110    4.157    4.157    5.048    5.046
 U,    kJ/mol =    2.541    1.239    1.239    1.779    1.770
 H,    kJ/mol =    2.541    1.239    1.239    1.779    1.770
 S,    J/K/mol=   13.428   17.763   17.562   17.293   17.093
 G,    kJ/mol =   -1.463   -4.057   -3.997   -3.377   -3.327
 U+ZPE,kJ/mol =    3.219    1.918    1.918    2.457    2.448
 H+ZPE,kJ/mol =    3.219    1.918    1.918    2.457    2.448
 G+ZPE,kJ/mol =   -0.785   -3.378   -3.318   -2.699   -2.648
 I(*)  - reduced moments of inertia are calculated by formula 1/I=1/Ia+1/Ib
 I(**) - reduced moments of inertia are calculated by Pitzer-Gwinn formulas
 Note: I(*)=I(**) for symmetric rotors.

Then, the total val­ues of ther­mo­dy­nam­ic func­tions of the mol­e­cule are print­ed togeth­er with total ener­gy val­ues cor­rect­ed for the total (electr+trans+rota+vibr+int.rot.) ther­mo­dy­nam­ic con­tri­bu­tions.

 Total values of thermodynamic functions (T=  298.15 K,  P= 101325.00 Pa) 
 taking into account  1 internal rotations  (in different approximations)
----------------------------------------------------------------------------------------------------------------------
                Harmonic   Free Rotator   Pitzer-Gwinn tables
               Ocsillator  I(*)    I(**)     I(*)    I(**)
----------------------------------------------------------------------------------------------------------------------
 Cv, J/K/mol=   77.176   73.223   73.223   74.114   74.111
 Cp, J/K/mol=   85.490   81.537   81.537   82.428   82.426
 U,  kJ/mol =   76.986   75.007   75.007   76.225   76.215
 H,  kJ/mol =   79.465   77.486   77.486   78.704   78.694
 S,  J/K/mol=  318.525  322.860  322.659  322.390  322.190
 G,  kJ/mol =  -15.503  -18.775  -18.715  -17.417  -17.367
----------------------------------------------------------------------------------------------------------------------
I(*)  - reduced moments of inertia are calculated by formula 1/I=1/Ia+1/Ib
I(**) - reduced moments of inertia are calculated by Pitzer-Gwinn formulas
Note: Internal rotation contributions to U, H, G include U0=H0=ZPE
in HO and PG models (U0=H0=0 in FR model)


 *** Total energy corrected by TD functions (including internal rotations) ****
 Minimum electronic-nuclear energy was taken from the input file at step:   4
 Minimum Etot:    -665.26660848  (Please check the program decision manually!)

          Harmonic             Free Rotator               Pitzer-Gwinn tables
         Ocsillator        I(*)           I(**)          I(*)           I(**)
E      -665.26660848  -665.26660848  -665.26660848  -665.26660848  -665.26660848
E+ZPE  -665.24329545  -665.24329545  -665.24329545  -665.24329545  -665.24329545
E+U    -665.23728531  -665.23803928  -665.23803928  -665.23757538  -665.23757887
E+H    -665.23634110  -665.23709507  -665.23709507  -665.23663117  -665.23663466
E+G    -665.27251350  -665.27375972  -665.27373686  -665.27324249  -665.27322325

The val­ues here are the results of total ener­gy cor­rect­ed by the con­tri­bu­tions of inter­nal rota­tion cal­cu­lat­ed at the cur­rent tem­per­a­ture and pres­sure ( pro­gram defaults: T=298.15K, P=101325Pa).

The file is fin­ished with the ther­mo­dy­nam­ic func­tions cal­cu­lat­ed at dif­fer­ent tem­per­a­tures. As before, the tem­per­a­ture range can be set man­u­al­ly (default 0–500 K). Now, the inter­nal rota­tion con­tri­bu­tions are includ­ed. The WARNING’s are just a result of numer­i­cal errors dur­ing the inter­po­la­tion of Pitzer-Gwinn tables for low tem­per­a­tures for light mol­e­cule (the tables are giv­en for the val­ues of Qf large enough and their extr­po­la­tion to zero is fre­quent­ly unro­bust). Just skip the val­ues with warn­ings. These val­ues can be incor­rect and should not be used. How­ev­er, they have no effect on the remain­ing val­ues.