The main goal of this work has been to obtain the Renner-Teller (RT) coupled-channel (CC) dynamics of the reaction, NH(a1Δ) + D(2S), considering the four following channels:

NH(a1Δ) + D(2S)--> N(2D) + HD(1Σg+) (depletion) (1)

NH(a1Δ) + D(2S) --> ND (a1Δ) + H(2S) (exchange) (2)

NH(a1Δ) + D(2S) --> NH(X3Σ-) + D(2S) (quenching) (3)

NH(a1Δ) + D(2S) --> ND(X3Σ-) + D(2S) (exchange+quenching) (4)

We have used the best available potential energy surfaces and have obtained initial-state-resolved probabilities, cross sections, rate constants and branching ratios via the time dependent real wave-packet method(1) (TDRWP) and flux analysis.

The two electronic states of NHD (X2B1 and A2A1) are the degenerate components of a linear 2Π state, thus giving rise to Renner-Teller(2) (RT) rovibronic nonadiabatic interactions allowing to change the electronic state. The RT effect is responsible to exchange and exchange + quenching reactions, so these reactions are not allowed under the Born-Oppenheimer approximation. Thus, the propagation of the RWP starts in the A2A1 excited potential energy surface(3) (PES), and when it reaches the H-N-D collinear arrangements it becomes possible to change the electronic state (jump into the X2B1 PES(3)) through the RT effect. Moreover, the NH2 X2B1 state has a deep minimum that traps the RWP for long time, increasing in this way the probability of nonadiabatic electronic transition (change of the electronic state).

Finally, the part of the RWP corresponding to reaction channels (1), (2), (3) and (4) is determined and the probability of each one of them is calculated, as a function of the initial conditions.

All the work detailed here is just the continuation of our previous project concerning the Renner-Teller dynamic and kinetic of atom+diatom reactions(5),(6).

References

(1) S. K. Gray and G. G. Balint-Kurti, J. Chem. Phys. 108, 950 (1998)

(2) C. Petrongolo, J. Chem. Phys. 89, 1297 (1988)

(3) Z.-W. Qu, H. Zhu, R. Schinke, L. Adam, and W. Hack, J. Chem. Phys. 122, 204313 (2005)

(4) S. Akpinar, P. Defazio, P. Gamallo and C. Petrongolo, J. Chem. Phys. 129, 174307 (2008)

(5) P. Gamallo and P. Defazio, J. Chem. Phys. 131, 44320 (2009)

(6) P. Defazio, P. Gamallo, M. González, S. Akpinar, B. Bussery-Honvault, P. Honvault and C. Petrongolo, J. Chem. Phys. 132, 104306-1 (2010)

The first thing we had to do was to test our parallel code in the CINECA machines. Some problems arise from this first step due to poor flexibility of IBM compiler included in the CINECA sp6 machine. Several test were performed and a lot of cpu time was wasted with the aim of solving the problem. Thus, we checked the code running some jobs and comparing the results with others obtained using other computers. The success of this step allowed us to begin the dynamic study of the reactions indicated in the previous section.We want to thank CINECA responsible to allow us to increase the cpu time for finishing all the work expected in the project.

Later on, we checked the convergence of the RT-CC-RWP (NH(a1Δ)+D(2S)) calculations verifying a large number of numerical parameters (e.g., rotational basis, mesh, number of iterations, etc.). Once the results were well converged, we performed the RWPs propagations for the title reaction.

The different initial conditions investigated have been the following: NH (v0=0, j0=2,3,4) and J=0,1,2,…,40, and K0=0,1,…,min(j,J). Because the CC method has been used, the RWP has been propagated so many times as given by the number of possible K final values (J+1), using a single processor for each propagation. A total of 80000 iteration steps are required to reach convergence and due to this the propagations of the RWPs have been very time demanding.

At present, we are carrying the analysis of all this amount of data using the flux and the asymptotic methods to obtain the probability of the four reaction channels. These probabilities will be the basis to obtain both the cross section and the rate constants for all processes. These rate constants will be compared with the experimental data available in the literature.

We simulate the interstellar medium (ISM) by solving numerically the 3-dimensional non-ideal magnetohydrodynamic (MHD) set of equations. The model includes the stellar gravity field in the solar neighbourhood, a stratified gas density, shearing due to differential rotation, radiative cooling and uv-heating.

For comparison we use simulations of purely hydrodynamic properties and then run them with the inclusion of a small seed magnetic field. The primary source of turbulence and heat is provided by supernovae (SNe). We include type I and type 2 SNe, which are exploded at rates similar to the observed rate in the solar neighbourhood and with a reasonably realistic randomised distribution by position.

The aim of this project is to produce a physically motivated dynamical structure for the ISM, which saturates to a quasi-steady state of turbulence, from which the typical temperature and density filling factors and turbulent velocities and vorticity can be calculated with respect to anisotropy and height of the ISM.

It is of interest to investigate the relationship between these variables and the rate of SNe. In addition the factors which determine the outward vertical velocity flows towards the galactic halo are of interest to understanding whether this constitutes a galactic wind or fountain.

In consideration of the magnetic field we investigate how these factors combine to influence the galactic dynamo. The mean galactic field appears to owe its structure to the rotational shear of the galaxy, whilst the random components appear to be strongly driven by the twisting and stretching from the persistent violent SNe events. Measurement and comparison of the mean and random components and there dependence on the rate of shear and SNe rate and distribution are of interest to understanding this process.

Within the ISM extreme variations within the density, temperature and velocity occur at close proximity. Combined with the violent shock fronts found in SNe remnants these properties are very challenging to model numerically, whilst retaining the essential physics effectively.

Before commencing the project, we had successfully modelled the essential properties described previously using type I SNe only and reaching simulation lengths of order 100-200 Myrs. As high temperatures or high velocities are generated by events in the simulation eventually the time step required to simulate accurately the processes becomes so small the simulation effectively stalls.

Our challenge with this project was to include the more numerous type II SNe and to extend the life span of the simulation to around a Gyr, whilst continuing to reproduce effectively and accurately a physically motivated galactic model. It can take 10-100 Myrs for turbulence to evolve to saturated steady state. To verify that results we obtain from our simulations are steady states and not numerical oscillations, we need to sustain these steady states over longer periods.

During the project we succeeded in solving many of the numerical challenges and have produced a simulation which exceeded 500Myrs and did not stall due to numerical difficulties, but because our budget for computing time has run out. This simulation included only hydrodynamics, although a similar run with magnetic field has also run over 100Myrs. In general it appears similar runs with magnetic field consistently can exceed those without, because the field appears to resist the propagation of excessive velocities or temperatures.

The simulations have produced results with interesting temperature and density filling factors and characteristic velocity patterns, which appear to be consistent with observations and current theory. To have a robust set of results, however, additional runs incorporating various alternative parameters and investigating the effect of higher resolution. We have already established that our current resolution of 5 parsecs cells is the minimum required to represent SNe expansion consistently, and this may yet need to be smaller to detect the dynamo properties consistently.

Complex long-term dynamics of reaction-diffusion problems requaires finer approach leading to quantitatively more reliable numerical schemes, for which the error estimate constant does not grow exponentially in time. One of approaches is known as the nonlinear Galerkin method and was introduced in M. Marion and R. Temam, Nonlinear Galerkin Methods, 1989. The method is rather general and allows several possibilities for the discretization of the problem in hand.

In our group we have been studying this method and its applications for the one dimensional problems for several year see e.g. J. Šembera and M. Beneš, Nonlinear Galerkin method for reaction-diffusion systems admitting invariant regions, 2001. In this paper, the method was applied to the numerical solution of the Brusselator reaction-diffusion systems and demonstrates its properties. Fourier spectral discretization and  Runge-Kutta method were used.

Later, we found the paper M. Marion and R. Temam, Nonlinar Galerkin Methods: The Finite Element Case, 1990, where the authors study application of the Nonlinear Galerkin approach for use with the finite element approximation in space. We applied this method to the numerical solution of selected reaction-diffcusion problems (the Gray-Scott model, the Brusselator model and the Phase-field model) in one spatial dimension, see e. g. J. Mach, Application of the Nonlinear Galerkin FEM to the solution of reaction-diffusion equations (under review), where we derive numerical scheme for the Gray-Scott reaction-diffusion system and compare this approach to the finite-difference method.

We continued the study and implemented the method for the solution of our test problem in two spatial dimensions. Some simplifications, which are possible in one spatial dimension, are no longer possible in higher dimensions. The computational time requirements of the serial implementation are high. For the numerical study we decided first to study the possibility of parallel implementation. The parallel implementation of the nonlinera Galerkin FEM method for selected reaction-diffusion system was the objetive. This will enable us to study the properties of the method in two spatial dimensions in a more convenient way.

During my stay at HLRS Stuttgart, Germany within the HPC-Europa2 program I studed the possibility of parallelization of our existing serial implementation of the nonlinear Galerkin FEM method for the solution of problem in two spatial dimensions.

The main part of my work here was devoted to implementation of selected parallel solver of linear systems of equations. We decided to try the MUMPS solver. It is a public domain software package for the multifrontal solution of large sparse linear systems on distributed memory computers. It can be used for a range of matrix types, so we will be able to use the experiance gained there even later, when working on different problem. We successfully implemented this solver into our code. For problem sizes we commonly use about 4x speedup was observed. We hope to gain more performace when having a closer look at the documentation later and use all this solver offers.