Anharmonicity plays a crucial role in determining the band profile and position of the interstellar Polycyclic Aromatic Hydrocarbons. Here we modeled the infrared spectrum of pyrene (C16H10) using \textrm{ab initio} simulations. We used our new AnharmoniCaOs code to describe the detailed structures of hot band transitions up to moderately high temperatures (600K), using Monte Carlo sampling of the states but beyond that due to high computational cost it is difficult to extend such fancy calculation up to very high temperature relevant for astrophysics. IR spectra of C16H10 beyond 600K was calculated using Density functional based tight binding molecular dynamics simulation. A gradual red shift of the band position and increasing band width is observed with increasing temperatures. The band positions at different temperatures were fitted with a linear regression and the anharmonicity factors were retrieved from the linear fits. Theoretical anharmonicity factors were compared to the recent laboratory results and previous gas phase analyses. This is an extension of a previous project, aimed at completing the calculations to obtain (relatively) high temperature anharmonic vibrational spectra of the pyrene (C16H10) Polycyclic Aromatic Hydrocarbon, using our AnharmoniCaOs code to describe the detailed structures of hot band transitions. We achieved our scientific goal, but the calculation was significantly more demanding than anticipated. The results have been submitted to the Journal of Molecular Spectroscopy by Chakraborty et al., I here describe the technical challenges we faced and the improvements we plan to overcome or mitigate them...

Chakraborty, S.

MULAS, Giacomo

OA@INAF - Istituto Nazionale di Astrofisica

2020Publication Year2022-03-29T10:56:26ZAcceptance in OA@INAFModelling anharmonic spectra of Polycyclic Aromatic Hydrocarbons at high temperaturesTitleMULAS, Giacomo; Chakraborty, S.Authorshttp://hdl.handle.net/20.500.12386/32014HandleMEMORIE DELLA SOCIETA ASTRONOMICA ITALIANAJournalMem. S.A.It. Vol. 91, 291c© SAIt 2020 Memorie dellaModelling anharmonic spectra of PolycyclicAromatic Hydrocarbons at high temperaturesG. Mulas1 and S. Chakraborty21 Istituto Nazionale di Astrofisica – Osservatorio Astronomico di Cagliari, Via della Scienza5, 09047 Selargius (CA), Italy, e-mail: giacomo.mulas@inaf.it2 Universitè de Toulouse – Institut de Recherche en Astrophysique et Planètologie, avenue duColonel Roche 9, 31028 Toulouse Cedex 4, FranceAbstract. Anharmonicity plays a crucial role in determining the band profile and positionof the interstellar Polycyclic Aromatic Hydrocarbons. Here we modeled the infrared spectrumof pyrene (C16H10) using ab initio simulations. We used our new AnharmoniCaOs code todescribe the detailed structures of hot band transitions up to moderately high temperatures(600K), using Monte Carlo sampling of the states but beyond that due to high computationalcost it is difficult to extend such fancy calculation up to very high temperature relevant forastrophysics. IR spectra of C16H10 beyond 600K was calculated using Density functional basedtight binding molecular dynamics simulation. A gradual red shift of the band position andincreasing band width is observed with increasing temperatures. The band positions at differenttemperatures were fitted with a linear regression and the anharmonicity factors were retrievedfrom the linear fits. Theoretical anharmonicity factors were compared to the recent laboratoryresults and previous gas phase analyses.This is an extension of a previous project, aimed at completing the calculations to obtain (rel-atively) high temperature anharmonic vibrational spectra of the pyrene (C16H10) PolycyclicAromatic Hydrocarbon, using our AnharmoniCaOs code to describe the detailed structuresof hot band transitions. We achieved our scientific goal, but the calculation was significantlymore demanding than anticipated. The results have been submitted to the Journal of MolecularSpectroscopy by Chakraborty et al., I here describe the technical challenges we faced and theimprovements we plan to overcome or mitigate them.1. IntroductionPolycyclic Aromatic Hydrocarbons (PAHs) area broad family of molecules composed offused aromatic rings, whose peripheral bondsare saturated by hydrogen atoms. PAHs (orclosely related species) are the most likely car-riers of the so-called Aromatic Infrared Bands(AIBs) at ∼3.3, 6.7, 7.7, 8.6, 11.3 µm (Leger1984; Allamandola 1985), which are amongthe strongest emission features observed inthe interstellar medium (Tielens 2008). PAHsin space are excited via electronic transitionsin the UV and visible, and re-emit most ofthis energy in the mid-IR in a vibrational de-excitation cascade. PAHs thus emit most of theflux in the AIBs from highly excited states cor-responding to vibrational temperatures of theorder of ∼1000 K (Draine 2011).Each AIB results from the superpositionof the emission of many different moleculesat many different temperatures, each of which292 Mulas: Modelling anharmonic spectra of PAHs at high Tinitial state(1, 1, 3, . . . , 0)initial state(1, 1, 3, . . . , 0)starting state(1, 0, 3, . . . , 0)starting state(2, 2, 3, . . . , 0)starting state(2, 1, 3, . . . , 0)transitionstransitioninitial state(1, 1, 3, . . . , 0)resonant state(2, 0, 3, . . . , 0)resonant state(0, 2, 3, . . . , 0)starting state(1, 0, 3, . . . , 0)resonant state(0, 1, 3, . . . , 0)resonant state(2, 1, 2, . . . , 0)starting state(2, 2, 3, . . . , 0)resonant state(1, 3, 3, . . . , 0)resonant state(3, 1, 3, . . . , 0)starting state(2, 1, 3, . . . , 0)resonant state(3, 0, 3, . . . , 0)resonant state(3, 3, 2, . . . , 0)resonancesresonancesresonancesresonancesinitial state(1, 1, 3, . . . , 0)resonant state(2, 0, 3, . . . , 0)resonant state(0, 2, 3, . . . , 0)starting state(1, 0, 3, . . . , 0)resonant state(0, 1, 3, . . . , 0)resonant state(2, 1, 2, . . . , 0)starting state(2, 2, 3, . . . , 0)resonant state(1, 3, 3, . . . , 0)resonant state(3, 1, 3, . . . , 0)starting state(2, 1, 3, . . . , 0)resonant state(3, 0, 3, . . . , 0)resonant state(3, 3, 2, . . . , 0)polyad 1polyad 3polyad 2Fig. 1. Schematic algorithm of polyad definitionbeing the superposition of many vibrationaltransitions from a statistical distribution of ex-cited states. Each of these individual transi-tions is slightly shifted by anharmonicity withrespect to the corresponding 1-0 fundamen-tal transition, resulting in complex, shifted,and broadened band profiles. Including an-harmonicity in AIB models is therefore cru-cial to link band positions and profiles to thephysics of the process of excitation and emis-sion. Our AnharmoniCaOs code 1 implementsthe generalized second-order vibrational per-turbation theory (GVPT2) (Gaw 1991; Martin1995; Barone 2005; Piccardo 2015) to com-pute anharmonic vibrational spectra from anyarbitrary vibrational state of a given molecule.We exploited this to obtain a Monte Carlo sam-pling of the vibrational states using the Wang-Landau approach, enabling us to obtain tem-perature dependent spectra of at finite, non-zero temperatures of Pyrene (C16H10) up to∼523 K.1 (https://sourceforge.net/projects/anharmonica/)This report covers both our previousproject CINECA-INAF B18 and the currentone, which is its continuation.2. Computational detailsWe made use of the Van Vleck approach toperturbation theory applied to molecular vi-brations (Sibert 1988a,b; McCoy 1992; Sibert2013). A simple recap can be found in ourpaper Mulas (2018). In our previous projectwe defined the details of the main approxi-mations and truncation strategies employed tomake the problem numerically tractable. Wealso validated method, truncation strategies,and thresholds against high-level calculationson smaller molecules. This has been publishedin our paper Mulas (2018), as well as cal-culations for Pyrene (C16H10) and Coronene(C24H12) at 0K, and comparisons with avail-able experimental data. In this report, we con-centrate on the use of our code on Pyrene to de-rive (relatively) high temperature spectra, morerelevant for astronomy.Mulas: Modelling anharmonic spectra of PAHs at high T 293Fig. 2. Left: density of harmonic vibrational states; right: density of anharmonic vibrational states.The core of the computational burden ofthis calculation can be divided in four maintasks:1. initialisation of the code, selection of res-onant Hamiltonian terms (to be used tobuild polyads and effective Hamiltonians innext step) and non-resonant terms, to be in-cluded in the perturbative treatment;2. construction of the polyads and of the cor-responding effective Hamiltonians;3. diagonalisation of the effectiveHamiltonians;4. calculation of the intensities of permittedtransitions between diagonalised states.Step 1 requires a relatively small amountof time, of the order of 1h on a single core.It could conceptually be parallelised in arelatively simple way, but it is not worththe effort, since it needs to be executed onlyonce and for all for a given molecule. It isimplemented as a separate code.Step 2, schematically described in Fig. 1,requires a relatively small amount of time, ofthe order of 1h on a single core when largepolyads are involved. It must be executed forevery step of the Wang-Landau random walkor, in general, for each state from which onewants to compute transitions. It is not paral-lelised, and it would not be trivial to do, sincethe list of states involved in all polyads wouldneed to be always available to all threads, thusrequiring some fairly complex locking andsynchronisations to avoid races. We found itmore convenient to also split this to a separatecode, which is run as a serial process.Step 3 is by far the most computationallyexpensive one. The effective Hamiltoniansare sparse matrices, and in principle onlya small subset of their eigensolutions areneeded, namely the ones which are necessaryto cover the vector subspace of the harmonicstates we are considering, in each specificiteration of the Wang-Landau random walk toobtain the Monte Carlo sampling of excitedstates. However, there are no available iterativeeigensolvers which are capable to selectivelyconverge on the eigensolutions we want: allknown algorithms for the iterative solutionof sparse eigenvalue problems select theeigensolutions on the basis of eigenvalues, notof eigenvectors. The implementation of theJacobi-Davidson method in the SLEPc library(Hernandez 2005, 2003) can be, to someextent, tweaked to favour the eigensolutionswe want, but still it turns out to be vastlyineffective compared to standard methods forthe complete solution of dense eigenvalueproblems. On the other hand, the latter workon dense matrices, hence they require hugeamounts of memory to accomodate them:a 100k × 100k matrix in double precisionrequires of the order of ∼100 Gbytes just tocontain it, the same amount for eigenvectors,and at least as much in working space, possiblyquite a bit more depending on the algorithm294 Mulas: Modelling anharmonic spectra of PAHs at high TFig. 3. Some of the strongest fundamental bands of Pyrene, computed for increasing temperatures. Black:14 K; red: 100 K; orange: 200 K; green: 300 K; blue: 423 K; magenta: 523 K.used, for a total of >300 Gbytes. This makesit mandatory to use parallel methods workingon distributed arrays, e. g. ScaLAPACK(Blackford 1997) or ELPA (Auckenthaler2011; Marek 2014). Both stock ScaLAPACKand ELPA have problems when dealing witharrays larger than ∼32k × 32k elements,and need to be compiled with long integerindices. An implementation of ScaLAPACKusing long integer indices is included in theIntel MKL library available on the CINECAmachines. We modified ourselves ELPA toinclude support for long integer indices. Thepdsyevx ScaLAPACK function turned out toalways return complete, correct results. ELPA,while in most cases about twice more efficientthan ScaLAPACK, sometimes unpredictablyfailed to converge gracefully issuing an errorcode, while a few times it just hung in endlessloops. We submitted error reports to theELPA developers, together with test codeto reproduce them. They acknowledged theproblems and replied that they would try to fixthem in forthcoming versions of their library.For the present runs, we therefore stuck toScaLAPACK’s pdsyevx function, employingup to 32 full nodes of galileo for the largestpolyads (mandated by RAM requirements).Step 4 is not computationally heavy, com-pared to the previous one. However, sinceeigenvectors are returned by the previous stepin distributed arrays, transition intensities arecomputed via calls to parallel BLACS func-tions, also implemented in the Intel MKL li-braries available on CINECA machines.Mulas: Modelling anharmonic spectra of PAHs at high T 295We should mention that, due to symmetryconsiderations, we partitioned Pyrene transi-tions in three separate groups that remain com-pletely separated, thereby reducing the size ofthe polyads and effective Hamiltonians.Summing up, therefore, the completeproject involved the following practical workon CINECA machines:1. running one serial code once for the initial-isation step 1 (repeated for each band type);2. running a serial code for each step inthe Wang-Landau random walk, to buildpolyads and effective Hamiltonians (step 2repeated for each band type, hence threetimes in total);3. running the parallel code for each stepin the Wang-Landau random walk, to di-agonalise the effective Hamiltonians andcompute transitions among resulting an-harmonic states (steps 3 and 4 repeatedfor each Wang-Landau step and each bandtype). Each of these parallel processes isrun on a number of nodes adequate forRAM requirements.The preliminary work of defining the stepsof the Wang-Landau random walk, as wellas post-processing to obtain energy-dependentand temperature-dependent spectra, were per-formed separately, on local workstations, sincethey do not require significant computationalresources and, conversely, do require signifi-cant interactive human intervention.We remark that determining RAM require-ments as a function of polyad size was not triv-ial: while the standard ScaLAPACK functionpdsyevx includes the possibility to issue a pre-liminary call that does no calculations but esti-mates the size of work space to be allocated,its implementation in the Intel MKL, for ef-ficiency reasons, dynamically allocates addi-tional memory in an undocumented amount,causing the code to just crash during the func-tion call (sometimes a couple of hours into it!)if the allocation fails. We thus had to prelim-inarly run a small subset of diagonalisations,and derive empirical rules as to the minimumamount of RAM (and hence of full nodes) thatis necessary to successfully complete the func-tion call for a given array size.3. High temperature calculations forPyreneUsing the Wang-Landau approach in the spaceof harmonic vibrational states, we obtained arelatively uniform sampling in (harmonic) en-ergy of the vibrational states of Pyrene upto about 11 × 103 cm−1. For each harmonicstate, we obtained a distribution of the an-harmonic states it contributes to, and theiranharmonic energies. Together with the ana-lytic density of harmonic states, this straight-forwardly yields the anharmonic density ofstates. Both harmonic and anharmonic densityof states are shown in Figure 2. As alreadyfound by previous works (Basire 2009; Parneix2013), the two densities are very similar, withthe anharmonic density curve running moreor less parallel to the harmonic one, slightlyhigher than it. We also obtained uniformlysampled energy-dependent vibrational spec-tra, from which we can obtain temperature-dependent spectra via a numerical Laplacetransform, weighted with the anharmonic den-sity of states. Since we sampled only ener-gies up to a given level, the numerical Laplacetransform is accurate only for T low enoughto give negligible weight to unsampled ener-gies, which in our case means that we can ob-tain reliable spectra only up to about ∼500 K.Figure 3 shows some of the strongest funda-mental bands of Pyrene, computed for a rangeof temperatures.The effect of temperature is apparent bothin terms of broadening and of position shift.The variation is also clearly stronger for somebands than for others. Some structure whichis easily resolved at 14 K is lost at the highesttemperatures, and plateaus, composed of verymany weak bands getting activated at hightemperatures, appear.Since AnharmoniCaOs also computesovertone, combination, and difference bands,they can also be quantitatively studied.Figure 4 shows two spectral regions in whichno fundamental bands are present.In the region between 2200 and 2400 cm−1no particular bands stand out, but a “carpet” ofweak features that, at high temperatures, tend296 Mulas: Modelling anharmonic spectra of PAHs at high TFig. 4. Combination bands of Pyrene, computed for increasing temperatures. Black: 14 K; red: 100 K;orange: 200 K; green: 300 K; blue: 423 K; magenta: 523 K.to blend in a pseudo-continuum. Conversely,between 5800 and 6000 cm−1 a group of bandsclearly stand out, corresponding to the energyof the C-H stretch overtones. Upon closer ex-amination of the calculation, these actuallyturn out to be dominated not by overtones ofthe IR-active C-H stretch modes (shown in thelowest panel of Figure 3, but of combinationbands of all the C-H stretch modes, includingIR-inactive ones. While their integrated totalintensity is relatively high, these bands appearto be very strongly broadened with increasingtemperatures. Since these bands would requirea PAH to be excitated to very high vibrationaltemperatures (T>1500 K) to achieve signifi-cant emission, if we extrapolate their behaviourthey are likely to be so shallow to be extremelydifficult to observe.4. Discussion and conclusionsWe showed that we can indeed quantitativelypredict the variation of vibrational band pro-files and positions with temperature. We arepreparing a publication in which we com-pare our theoretical calculation with avail-able experimental data, to assess their accu-racy (Chakraborty et al. 2021). It is clear frompreliminary results that we predict the correctqualitative behaviour, but we still need to checkhow precisely the calculations match the ex-periment. As mentioned in the report for theprevious project, we are hindered by the lackof eigensolvers tailored to our problem. Weare however in the process of rewriting ourAnharmoniCaOs code to make it more effi-cient in two main ways. The first one is toMulas: Modelling anharmonic spectra of PAHs at high T 297add a step of sorting, merging, and reorder-ing of the polyads before building the effectiveHamiltonians. In this way, if a given polyadappears several times, e. g. once as the lowerstate of a transition and another as the up-per state of another transition, it will be com-puted only once. The second planned improve-ment is to try to take advantage of the factthat the states, and the polyads, change grad-ually along the Wang-Landau random walkamong harmonic states. As a consequence,also the corresponding effective Hamiltoniansare relatively close between consecutive stepsof the random walk. One may therefore usethe solutions of the previous step as a “start-ing point” for the next one, and use an iter-ative solver instead of the standard RelativelyRobust Representations (Dhillon 1997) imple-mented e. g. in ScaLAPACK. It remains to besee how effective these changes will prove tobe from the computational point of view.Acknowledgements. We acknowledge the comput-ing centre of Cineca and INAF, under the coordi-nation of the “Accordo Quadro MoU per lo svolgi-mento di attività congiunta di ricerca Nuove fron-tiere in Astrofisica: HPC e Data Exploration dinuova generazione”, for the availability of comput-ing resources and support.ReferencesAllamandola, L. J., Tielens, A. G. G. M.,Barker, J. R. 1985, ApJ, 290, L25Auckenthaler, T., et al. 2011, Parallel Comput.,37, 783Barone, V. 2005, J. Chem. Phys., 122, 014108Basire, M., et al. 2009, J. Phys. Chem. A, 113,6947Blackford, et al. 1997, ScaLAPACK Users’Guide http://www.netlib.org/scalapack/slug/Chakraborty, S., Mulas, G., Rapacioli, M., etal. 2021, arXiv:2102.06582Dhillon, I. 1997, Report n. UCB/CSD-97-971,EECS Department, University of California,BerkeleyDraine, B.T. 2011, EAS Publications Series,46, 29Gaw, F., et al. 1991, in Advances in MolecularVibrations and Collision Dynamics: AResearch Annual, ed. J. M. Bowman and M.A. Ratner (JAI Press, Greenwich, Conn.),1B, 169Hernandez, V., Vidal, V. 2003, Lect. NotesComput. Sci., 2565, 377Hernandez, V., Roman, J. E., Vidal, V. 2005,ACM Trans. Math. Software, 31, 351Léger, A., Puget, J. L. 1984, A&A, 137, L5Marek, A., et al. 2014, J. Phys. - Condens.Mat., 26, 213201Martin, J. M. L., et al. 1995, J. Chem. Phys.,103, 2589McCoy, A. B., Sibert, E. L. 1992, Mol. Phys.,77, 697Mulas, G., et al. 2018, J. Chem. Phys., 149,144102Parneix, P., Basire, M., Calvo, F. 2013, J. Phys.Chem. A, 117, 3954Piccardo, M., Bloino, J., Barone, V. 2015, Int.J. Quantum Chem., 115, 948Sibert, E. L. 1988a, J. Chem. Phys., 88, 4378Sibert, E. L. 1988b, Comput. Phys. Commun.,51, 149Sibert, E. L. 2013, Mol. Phys., 111, 2093Tielens, A. G. G. M. 2008, ARA&A, 46, 289

Modelling anharmonic spectra of Polycyclic Aromatic Hydrocarbons at high temperatures

https://openaccess.inaf.it/bitstream/20.500.12386/32014/1/291.pdf

Modelling anharmonic spectra of Polycyclic Aromatic Hydrocarbons at high temperatures

Abstract

Similar works

Full text

Available Versions

OA@INAF - Istituto Nazionale di Astrofisica