`arkane` unable to infer molecular weights with pdep

[liqiang4567]

Topic

General area which your question is related to.

• [ ] Arkane (formerly CanTherm)

Question

I want to obtain the pressure-dependent rate coefficients of C2H2 + C5H5 = C7H7 through Arkane and obtain the relevant information through the gaussian software. An IndexError: list index out of range error occurred

Automated Reaction Kinetics and Network Exploration (Arkane) Version: 3.2.0 RMG Installed via docker

input.py as follows:

title = 'Example calculation of C2H2 + C5H5 = C7H7'

description = \
"""
This example is for a rigid rotor harmonic oscillator TST calculation calclated at the 
UWB97XD/def2TZVP level of DFT with a tight grid.  Tunneling is not included in this calculation.
"""
#modelChemistry = LevelOfTheory(method='CCSD(T)-F12', basis='cc-pVTZ-F12', software='molpro')
#modelChemistry = "WB97XD/def2TZVP"

modelChemistry = LevelOfTheory(method="wb97xd", basis="augccpvtz",software="gaussian")

frequencyScaleFactor = 0.955
useHinderedRotors = False
useBondCorrections = False

species('C2H2', 'C#C.py')
species('C5H5', 'C1=C[CH]C=C1.py')

transitionState('ts1', 'ts1.py')
#species('C7H7', '[CH]=CC1C=CC=C1.py',collisionModel = TransportData(sigma=(6.03,'angstrom'), epsilon=(3667.12,'J/mol')))
species('C7H7', '[CH]=CC1C=CC=C1.py')

species(
    label = 'Ar',
    structure = SMILES('[Ar]'),
    E0 = (-6.19426,'kJ/mol'),
    spinMultiplicity = 1,
    opticalIsomers = 1,
    molecularWeight = (39.348,'amu'),
    collisionModel = TransportData(shapeIndex=0, epsilon=(1134.93,'J/mol'), sigma=(3.33,'angstroms'), dipoleMoment=(0,'C*m'), polarizability=(0,'angstroms^3'), rotrelaxcollnum=0.0, comment="""GRI-Mech"""),
    energyTransferModel = SingleExponentialDown(alpha0=(133,'cm^-1'), T0=(300,'K'), n=0.85),
    thermo = NASA(polynomials=[NASAPolynomial(coeffs=[2.5,2.64613e-14,-3.72536e-17,1.7192e-20,-2.44483e-24,-745,4.3663], Tmin=(100,'K'), Tmax=(3802.52,'K')), NASAPolynomial(coeffs=[2.5,1.04239e-10,-3.81845e-14,6.18592e-18,-3.73869e-22,-745,4.3663], Tmin=(3802.52,'K'), Tmax=(5000,'K'))], Tmin=(100,'K'), Tmax=(5000,'K'), E0=(-6.19426,'kJ/mol'), Cp0=(20.7862,'J/(mol*K)'), CpInf=(20.7862,'J/(mol*K)'), label="""Ar""", comment="""Thermo library: BurkeH2O2"""),
)

thermo('C2H2','NASA')
thermo('C5H5','NASA')
thermo('C7H7','NASA')

reaction(
    label = 'C2H2 + C5H5 <=> C7H7',
    reactants = ['C2H2', 'C5H5'],
    products = ['C7H7'], #product channels are only important if tunneling is computed
    transitionState = 'ts1',
    tunneling='Eckart', #we dont want to comput tunneling for this calculation
)

network(
    label = 'C2H2_net + C5H5_net <=> C7H7_net',
    isomers = [
        'C7H7',
    ],
    bathGas = {
        'Ar': 1.0,
    }
)

################################################################################

pressureDependence(
    'C2H2_net + C5H5_net <=> C7H7_net',
    #Tlist = ([300, 400, 600, 800, 1000, 1250, 1500, 1750, 2000],'K')
    #Plist = ([0.01, 0.1, 1.0, 10.0, 100.0],'bar')
    Tlist = ([500, 600,700, 800, 900,950, 1000],'K'),
    Plist = ([10.0,50.0,70.0,100.0],'atm'),
    #Tmin=(300.0,'K'), Tmax=(2000.0,'K'), Tcount=8,
    #Pmin=(0.01,'bar'), Pmax=(100.0,'bar'), Pcount=5,
    maximumGrainSize = (0.5,'kcal/mol'),
    minimumGrainCount = 200,
    method = 'modified strong collision',
    interpolationModel = ('chebyshev', 6, 4),
    #activeKRotor = True, 
    activeJRotor = True,
)

Output is as follows:
################################################################
#                                                              #
# Automated Reaction Kinetics and Network Exploration (Arkane) #
#                                                              #
#   Version: 3.2.0                                             #
#   Authors: RMG Developers (rmg_dev@mit.edu)                  #
#   P.I.s:   William H. Green (whgreen@mit.edu)                #
#            Richard H. West (r.west@neu.edu)                  #
#   Website: http://reactionmechanismgenerator.github.io/      #
#                                                              #
################################################################

The current git HEAD for RMG-Py is:
    036ab3f8ca0f94f567b50b0b83110bab0a14a35f
    Wed Aug 16 15:56:48 2023 -0400

The current git HEAD for RMG-database is:
    b7ff16364a07c9a51a34303aa28407a83455a3e4
    Tue Aug 8 09:37:08 2023 -0400

Databases not found. Making databases
Loading frequencies library from halogens_G4.py in /rmg/RMG-database/input/statmech/libraries...
Loading frequencies group database from /rmg/RMG-database/input/statmech/groups...
Loading transport library from GRI-Mech.py in /rmg/RMG-database/input/transport/libraries...
Loading transport library from NIST_Fluorine.py in /rmg/RMG-database/input/transport/libraries...
Loading transport library from NOx2018.py in /rmg/RMG-database/input/transport/libraries...
Loading transport library from OneDMinN2.py in /rmg/RMG-database/input/transport/libraries...
Loading transport library from PrimaryTransportLibrary.py in /rmg/RMG-database/input/transport/libraries...
Loading transport group database from /rmg/RMG-database/input/transport/groups...
Loading species C2H2...
Loading species C5H5...
Loading transition state ts1...
Loading species C7H7...
Loading species Ar...
Loading reaction C2H2 + C5H5 <=> C7H7...
Loading network C2H2_net + C5H5_net <=> C7H7_net...

Loading statistical mechanics parameters for C2H2...
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: Dinfh
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: Dinfh
Saving statistical mechanics parameters for C2H2...
Loading statistical mechanics parameters for C5H5...
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: C2v
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: C2v
Saving statistical mechanics parameters for C5H5...
Loading statistical mechanics parameters for ts1...
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: Cs
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: Cs
Saving statistical mechanics parameters for ts1...
Loading statistical mechanics parameters for C7H7...
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: Cs
Symmetry input file written to /rmg/RMG-Py/C5H5_C2H2/Rnx28/scratch/0.symm
Point group: Cs
Saving statistical mechanics parameters for C7H7...
Saving thermo for C2H2...
Saving thermo for C5H5...
Saving thermo for C7H7...
========================================================================
C2H2_net + C5H5_net <=> C7H7_net network information
----------------------------------------------------
Isomers:
    C7H7                                                  450.944 kJ/mol
Reactant channels:
Product channels:
    C2H2 + C5H5                                           474.797 kJ/mol
Path reactions:
    C2H2 + C5H5 <=> C7H7                                  539.186 kJ/mol
Net reactions:
========================================================================

Calculating densities of states for C2H2_net + C5H5_net <=> C7H7_net network...
Using 332 grains from 450.94 to 874.31 kJ/mol in steps of 1.28 kJ/mol to compute densities of states
Calculating phenomenological rate coefficients for C2H2 + C5H5 <=> C7H7...
Using 200 grains from 450.94 to 705.48 kJ/mol in steps of 1.28 kJ/mol to compute the k(T,P) values at 500 K
Traceback (most recent call last):
  File "../../Arkane.py", line 59, in <module>
    arkane.execute()
  File "/rmg/RMG-Py/arkane/main.py", line 284, in execute
    job.execute(output_file=output_file, plot=self.plot)
  File "/rmg/RMG-Py/arkane/pdep.py", line 290, in execute
    self.K = self.network.calculate_rate_coefficients(self.Tlist.value_si, self.Plist.value_si, self.method)
  File "/rmg/RMG-Py/rmgpy/pdep/network.py", line 268, in calculate_rate_coefficients
    self.set_conditions(T, P)
  File "/rmg/RMG-Py/rmgpy/pdep/network.py", line 449, in set_conditions
    self.calculate_collision_model()
  File "/rmg/RMG-Py/rmgpy/pdep/network.py", line 907, in calculate_collision_model
    coll_freq[i] = isomer.calculate_collision_frequency(self.T, self.P, self.bath_gas)
  File "rmgpy/pdep/configuration.pyx", line 174, in rmgpy.pdep.configuration.Configuration.calculate_collision_frequency
  File "rmgpy/pdep/configuration.pyx", line 188, in rmgpy.pdep.configuration.Configuration.calculate_collision_frequency
  File "rmgpy/species.py", line 810, in rmgpy.species.Species.get_transport_data
  File "rmgpy/species.py", line 796, in rmgpy.species.Species.generate_transport_data
IndexError: list index out of range

[JacksonBurns]

@xiaoruiDong have you seen something like this before?

[xiaoruiDong]

I haven't. The error is raised on this line.

https://github.com/ReactionMechanismGenerator/RMG-Py/blob/036ab3f8ca0f94f567b50b0b83110bab0a14a35f/rmgpy/species.py#L796C9-L796C51 if isinstance(self.molecule[0], Molecule):

This implies that the species object has no defined corresponding molecular structure. Without such information, RMG cannot estimate its transport properties. I think the problem can be solved by providing a molecular identifier (either SMILES or InChI) when defining the species in the input file. For example:

species('C2H2', 'C#C.py', structure=SMILES('C#C'))
species('C5H5', 'C1=C[CH]C=C1.py', structure=SMILES('C1=C[CH]C=C1'))
species('C7H7', '[CH]=CC1C=CC=C1.py', structure=SMILES('[CH]=CC1C=CC=C1'))

@liqiang4567 Can you try this out and report if it solves your problem? Thanks.

[liqiang4567]

Thank you very much for your answer, the program is running normally. But I have one doubt:

1. Why do I need structure or transport data when calculating pressure-dependent rate coefficients?
2. How does this affect my results? I want to get the speed all through quantum mechanics。

[xiaoruiDong]

Transport data (such as Lennard-Jones parameters) are needed to compute the collision frequency. RMG uses molecular structure to estimate the transport data using Joback group additivity (see RMG database paper for more details). You can compute them instead of estimating them, e.g., through the approach elaborated in this paper, but it requires extra calculations beyond the QM calculation you have done.

[liqiang4567]

Thank you for your reply. If I can get the Transport data, I don't need the structure. Can you write it like this?

species('C2H2', 'C#C.py',collisionModel = TransportData(sigma=(4.35,'angstrom'), epsilon=(2353.39,'J/mol')))
species('C5H5', 'C1=C[CH]C=C1.py',collisionModel = TransportData(sigma=(5.71,'angstrom'), epsilon=(3411.73,'J/mol')))
species('C7H7', '[CH]=CC1C=CC=C1.py',collisionModel = TransportData(sigma=(6.03,'angstrom'), epsilon=(3667.12,'J/mol')))

But there seems to be a mistake：

Calculating densities of states for C2H2_net + C5H5_net <=> C7H7_net network...
Using 1358 grains from 450.94 to 871.80 kJ/mol in steps of 0.31 kJ/mol to compute densities of states
Calculating phenomenological rate coefficients for C2H2 + C5H5 <=> C7H7...
Using 500 grains from 450.94 to 605.70 kJ/mol in steps of 0.31 kJ/mol to compute the k(T,P) values at 200 K
Traceback (most recent call last):
  File "../../Arkane.py", line 59, in <module>
    arkane.execute()
  File "/rmg/RMG-Py/arkane/main.py", line 284, in execute
    job.execute(output_file=output_file, plot=self.plot)
  File "/rmg/RMG-Py/arkane/pdep.py", line 290, in execute
    self.K = self.network.calculate_rate_coefficients(self.Tlist.value_si, self.Plist.value_si, self.method)
  File "/rmg/RMG-Py/rmgpy/pdep/network.py", line 268, in calculate_rate_coefficients
    self.set_conditions(T, P)
  File "/rmg/RMG-Py/rmgpy/pdep/network.py", line 449, in set_conditions
    self.calculate_collision_model()
  File "/rmg/RMG-Py/rmgpy/pdep/network.py", line 907, in calculate_collision_model
    coll_freq[i] = isomer.calculate_collision_frequency(self.T, self.P, self.bath_gas)
  File "rmgpy/pdep/configuration.pyx", line 174, in rmgpy.pdep.configuration.Configuration.calculate_collision_frequency
  File "rmgpy/pdep/configuration.pyx", line 200, in rmgpy.pdep.configuration.Configuration.calculate_collision_frequency
AttributeError: 'NoneType' object has no attribute 'value_si'

[xiaoruiDong]

From what I can tell, it complains about missing molecular weights https://github.com/ReactionMechanismGenerator/RMG-Py/blob/036ab3f8ca0f94f567b50b0b83110bab0a14a35f/rmgpy/pdep/configuration.pyx#L200

I think this is an unexpected behavior, as it should be inferred from the calculation. Anyway, can you try explicitly assigning molecular weights to the species as a temporary workaround?

[liqiang4567]

I have added molecularWeight, and even though this should be able to be obtained from the calculation, an error still appears:

Calculating densities of states for C2H2_net + C5H5_net <=> C7H7_net network...
Using 1358 grains from 450.94 to 871.80 kJ/mol in steps of 0.31 kJ/mol to compute densities of states
Calculating phenomenological rate coefficients for C2H2 + C5H5 <=> C7H7...
Using 500 grains from 450.94 to 605.70 kJ/mol in steps of 0.31 kJ/mol to compute the k(T,P) values at 200 K
Using 500 grains from 450.94 to 672.22 kJ/mol in steps of 0.44 kJ/mol to compute the k(T,P) values at 400 K
Using 500 grains from 450.94 to 705.48 kJ/mol in steps of 0.51 kJ/mol to compute the k(T,P) values at 500 K
Using 500 grains from 450.94 to 738.73 kJ/mol in steps of 0.58 kJ/mol to compute the k(T,P) values at 600 K
Using 500 grains from 450.94 to 771.99 kJ/mol in steps of 0.64 kJ/mol to compute the k(T,P) values at 700 K
Using 500 grains from 450.94 to 805.25 kJ/mol in steps of 0.71 kJ/mol to compute the k(T,P) values at 800 K
Using 500 grains from 450.94 to 838.51 kJ/mol in steps of 0.78 kJ/mol to compute the k(T,P) values at 900 K
Using 500 grains from 450.94 to 871.76 kJ/mol in steps of 0.84 kJ/mol to compute the k(T,P) values at 1000 K
Saving pressure dependence results for network C2H2_net + C5H5_net <=> C7H7_net...
Traceback (most recent call last):
  File "../../Arkane.py", line 59, in <module>
    arkane.execute()
  File "/rmg/RMG-Py/arkane/main.py", line 284, in execute
    job.execute(output_file=output_file, plot=self.plot)
  File "/rmg/RMG-Py/arkane/pdep.py", line 295, in execute
    self.save(output_file)
  File "/rmg/RMG-Py/arkane/pdep.py", line 473, in save
    f.write("#" + spc.label + "  SMILES: " + spc.molecule[0].to_smiles() + "\n")
IndexError: list index out of range

[xiaoruiDong]

Then I would suggest adding structure information back (e.g., structure=SMILES('C#C')), as indicated by documentation (structure is required) and the error. It is needed to write the output files. Since you have provided the collisionModel, Arkane will only use the structure information to write output files.

Thank you for catching the issues. We will fix the bug related to molecular weights and update the documentation accordingly.

[liqiang4567]

Ok, thank you very much for your reply

[JacksonBurns]

@xiaoruiDong which species were actually causing the problem (Arkane not being able to get their Molecular Weight on its own)?

[github-actions[bot]]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.