Novel algorithms for electronic structure based molecular dynamics with linear system-size scaling

Die vorliegende Arbeit behandelt die Entwicklung und Verbesserung von linear skalierenden Algorithmen für Elektronenstruktur basierte Molekulardynamik. Molekulardynamik ist eine Methode zur Computersimulation des komplexen Zusammenspiels zwischen Atomen und Molekülen bei endlicher Temperatur. Ein entscheidender Vorteil dieser Methode ist ihre hohe Genauigkeit und Vorhersagekraft. Allerdings verhindert der Rechenaufwand, welcher grundsätzlich kubisch mit der Anzahl der Atome skaliert, die Anwendung auf große Systeme und lange Zeitskalen. Ausgehend von einem neuen Formalismus, basierend auf dem großkanonischen Potential und einer Faktorisierung der Dichtematrix, wird die Diagonalisierung der entsprechenden Hamiltonmatrix vermieden. Dieser nutzt aus, dass die Hamilton- und die Dichtematrix aufgrund von Lokalisierung dünn besetzt sind. Das reduziert den Rechenaufwand so, dass er linear mit der Systemgröße skaliert. Um seine Effizienz zu demonstrieren, wird der daraus entstehende Algorithmus auf ein System mit flüssigem Methan angewandt, das extremem Druck (etwa 100 GPa) und extremer Temperatur (2000 - 8000 K) ausgesetzt ist. In der Simulation dissoziiert Methan bei Temperaturen oberhalb von 4000 K. Die Bildung von sp²-gebundenem polymerischen Kohlenstoff wird beobachtet. Die Simulationen liefern keinen Hinweis auf die Entstehung von Diamant und wirken sich daher auf die bisherigen Planetenmodelle von Neptun und Uranus aus. Da das Umgehen der Diagonalisierung der Hamiltonmatrix die Inversion von Matrizen mit sich bringt, wird zusätzlich das Problem behandelt, eine (inverse) p-te Wurzel einer gegebenen Matrix zu berechnen. Dies resultiert in einer neuen Formel für symmetrisch positiv definite Matrizen. Sie verallgemeinert die Newton-Schulz Iteration, Altmans Formel für beschränkte und nicht singuläre Operatoren und Newtons Methode zur Berechnung von Nullstellen von Funktionen. Der Nachweis wird erbracht, dass die Konvergenzordnung immer mindestens quadratisch ist und adaptives Anpassen eines Parameters q in allen Fällen zu besseren Ergebnissen führt.

Molecular dynamics simulations of silicate and borate glasses and melts : structure, diffusion dynamics and vibrational properties

Molecular dynamics simulations of silicate and borate glasses and melts: Structure, diffusion dynamics and vibrational properties. In this work computer simulations of the model glass formers SiO2 and B2O3 are presented, using the techniques of classical molecular dynamics (MD) simulations and quantum mechanical calculations, based on density functional theory (DFT). The latter limits the system size to about 100−200 atoms. SiO2 and B2O3 are the two most important network formers for industrial applications of oxide glasses. Glass samples are generated by means of a quench from the melt with classical MD simulations and a subsequent structural relaxation with DFT forces. In addition, full ab initio quenches are carried out with a significantly faster cooling rate. In principle, the structural properties are in good agreement with experimental results from neutron and X-ray scattering, in all cases. A special focus is on the study of vibrational properties, as they give access to low-temperature thermodynamic properties. The vibrational spectra are calculated by the so-called ”frozen phonon” method. In all cases, the DFT curves show an acceptable agreement with experimental results of inelastic neutron scattering. In case of the model glass former B2O3, a new classical interaction potential is parametrized, based on the liquid trajectory of an ab initio MD simulation at 2300 K. In this course, a structural fitting routine is used. The inclusion of 3-body angular interactions leads to a significantly improved agreement of the liquid properties of the classical MD and ab initio MD simulations. However, the generated glass structures, in all cases, show a significantly lower fraction of 3-membered planar boroxol rings as predicted by experimental results (f=60%-80%). The largest boroxol ring fraction of f=15±5% is observed in the full ab initio quenches from 2300 K. In case of SiO2, the glass structures after the quantum mechanical relaxation are the basis for calculations of the linear thermal expansion coefficient αL(T), employing the quasi-harmonic approximation. The striking observation is a change change of sign of αL(T) going along with a temperature range of negative αL(T) at low temperatures, which is in good agreement with experimental results.

Testing of subgrid scale (SGS) models for large-eddy simulation (LES) of turbulent channel flow

Sub-grid scale (SGS) models are required in order to model the influence of the unresolved small scales on the resolved scales in large-eddy simulations (LES), the flow at the smallest scales of turbulence. In the following work two SGS models are presented and deeply analyzed in terms of accuracy through several LESs with different spatial resolutions, i.e. grid spacings. The first part of this thesis focuses on the basic theory of turbulence, the governing equations of fluid dynamics and their adaptation to LES. Furthermore, two important SGS models are presented: one is the Dynamic eddy-viscosity model (DEVM), developed by \cite{germano1991dynamic}, while the other is the Explicit Algebraic SGS model (EASSM), by \cite{marstorp2009explicit}. In addition, some details about the implementation of the EASSM in a Pseudo-Spectral Navier-Stokes code \cite{chevalier2007simson} are presented. The performance of the two aforementioned models will be investigated in the following chapters, by means of LES of a channel flow, with friction Reynolds numbers $Re_\tau=590$ up to $Re_\tau=5200$, with relatively coarse resolutions. Data from each simulation will be compared to baseline DNS data. Results have shown that, in contrast to the DEVM, the EASSM has promising potentials for flow predictions at high friction Reynolds numbers: the higher the friction Reynolds number is the better the EASSM will behave and the worse the performances of the DEVM will be. The better performance of the EASSM is contributed to the ability to capture flow anisotropy at the small scales through a correct formulation for the SGS stresses. Moreover, a considerable reduction in the required computational resources can be achieved using the EASSM compared to DEVM. Therefore, the EASSM combines accuracy and computational efficiency, implying that it has a clear potential for industrial CFD usage.

Impact of irrigation in the Po valley (Italy) on meteorological parameters through numerical simulation

Recent studies found that soil-atmosphere coupling features, through soil moisture, have been crucial to simulate well heat waves amplitude, duration and intensity. Moreover, it was found that soil moisture depletion both in Winter and Spring anticipates strong heat waves during the Summer. Irrigation in geophysical studies can be intended as an anthropogenic forcing to the soil-moisture, besides changes in land proprieties. In this study, the irrigation was add to a LAM hydrostatic model (BOLAM) and coupled with the soil. The response of the model to irrigation perturbation is analyzed during a dry Summer season. To identify a dry Summer, with overall positive temperature anomalies, an extensive climatological characterization of 2015 was done. The method included a statistical validation on the reference period distribution used to calculate the anomalies. Drought conditions were observed during Summer 2015 and previous seasons, both on the analyzed region and the Alps. Moreover July was characterized as an extreme event for the referred distribution. The numerical simulation consisted on the summer season of 2015 and two run: a control run (CTR), with the soil coupling and a perturbed run (IPR). The perturbation consists on a mask of land use created from the Cropland FAO dataset, where an irrigation water flux of 3 mm/day was applied from 6 A.M. to 9 A.M. every day. The results show that differences between CTR and IPR has a strong daily cycle. The main modifications are on the air masses proprieties, not on to the dynamics. However, changes in the circulation at the boundaries of the Po Valley are observed, and a diagnostic spatial correlation of variable differences shows that soil moisture perturbation explains well the variation observed in the 2 meters height temperature and in the latent heat fluxes.On the other hand, does not explain the spatial shift up and downslope observed during different periods of the day. Given the results, irrigation process affects the atmospheric proprieties on a larger scale than the irrigation, therefore it is important in daily forecast, particularly during hot and dry periods.

Pi release from myosin: a simulation analysis of possible pathways

The release of phosphate (Pi) is an important element in actomyosin function and has been shown to be accelerated by the binding of myosin to actin. To provide information about the structural elements important for Pi release, possible escape pathways from various isolated myosin II structures have been determined by molecular dynamics simulations designed for studying such slow processes. The residues forming the pathways were identified and their role was evaluated by mutant simulations. Pi release is slow in the pre-powerstroke structure, an important element in preventing the powerstroke prior to actin binding, and is much more rapid for Pi modeled into the post-rigor and rigor-like structures. The previously proposed backdoor route is dominant in the pre-powerstroke and post-rigor states, whereas a different path is most important in the rigor-like state. This finding suggests a mechanism for the actin-activated acceleration of Pi release.

Simulating the impact of four control strategies on the population dynamics of Neospora caninum infection in Swiss dairy cattle

A dynamic deterministic simulation model was developed to assess the impact of different putative control strategies on the seroprevalence of Neospora caninum in female Swiss dairy cattle. The model structure comprised compartments of "susceptible" and "infected" animals (SI-model) and the cattle population was divided into 12 age classes. A reference model (Model 1) was developed to simulate the current (status quo) situation (present seroprevalence in Switzerland 12%), taking into account available demographic and seroprevalence data of Switzerland. Model 1 was modified to represent four putative control strategies: testing and culling of seropositive animals (Model 2), discontinued breeding with offspring from seropositive cows (Model 3), chemotherapeutic treatment of calves from seropositive cows (Model 4), and vaccination of susceptible and infected animals (Model 5). Models 2-4 considered different sub-scenarios with regard to the frequency of diagnostic testing. Multivariable Monte Carlo sensitivity analysis was used to assess the impact of uncertainty in input parameters. A policy of annual testing and culling of all seropositive cattle in the population reduced the seroprevalence effectively and rapidly from 12% to <1% in the first year of simulation. The control strategies with discontinued breeding with offspring from all seropositive cows, chemotherapy of calves and vaccination of all cattle reduced the prevalence more slowly than culling but were still very effective (reduction of prevalence below 2% within 11, 23 and 3 years of simulation, respectively). However, sensitivity analyses revealed that the effectiveness of these strategies depended strongly on the quality of the input parameters used, such as the horizontal and vertical transmission factors, the sensitivity of the diagnostic test and the efficacy of medication and vaccination. Finally, all models confirmed that it was not possible to completely eradicate N. caninum as long as the horizontal transmission process was not interrupted.

Reynolds-stress turbulence model in the KIVA code for engine simulation

A Reynolds-Stress Turbulence Model has been incorporated with success into the KIVA code, a computational fluid dynamics hydrocode for three-dimensional simulation of fluid flow in engines. The newly implemented Reynolds-stress turbulence model greatly improves the robustness of KIVA, which in its original version has only eddy-viscosity turbulence models. Validation of the Reynolds-stress turbulence model is accomplished by conducting pipe-flow and channel-flow simulations, and comparing the computed results with experimental and direct numerical simulation data. Flows in engines of various geometry and operating conditions are calculated using the model, to study the complex flow fields as well as confirm the model’s validity. Results show that the Reynolds-stress turbulence model is able to resolve flow details such as swirl and recirculation bubbles. The model is proven to be an appropriate choice for engine simulations, with consistency and robustness, while requiring relatively low computational effort.

COMPUTATIONAL ANALYSIS OF THE AEROELASTIC DYNAMICS OF A WIND TURBINE BLADE WITH VARIABLE-SPEED STALL CONTROL

The accuracy of simulating the aerodynamics and structural properties of the blades is crucial in the wind-turbine technology. Hence the models used to implement these features need to be very precise and their level of detailing needs to be high. With the variety of blade designs being developed the models should be versatile enough to adapt to the changes required by every design. We are going to implement a combination of numerical models which are associated with the structural and the aerodynamic part of the simulation using the computational power of a parallel HPC cluster. The structural part models the heterogeneous internal structure of the beam based on a novel implementation of the Generalized Timoshenko Beam Model Technique.. Using this technique the 3-D structure of the blade is reduced into a 1-D beam which is asymptotically equivalent. This reduces the computational cost of the model without compromising its accuracy. This structural model interacts with the Flow model which is a modified version of the Blade Element Momentum Theory. The modified version of the BEM accounts for the large deflections of the blade and also considers the pre-defined structure of the blade. The coning, sweeping of the blade, tilt of the nacelle and the twist of the sections along the blade length are all computed by the model which aren’t considered in the classical BEM theory. Each of these two models provides feedback to the other and the interactive computations lead to more accurate outputs. We successfully implemented the computational models to analyze and simulate the structural and aerodynamic aspects of the blades. The interactive nature of these models and their ability to recompute data using the feedback from each other makes this code more efficient than the commercial codes available. In this thesis we start off with the verification of these models by testing it on the well-known benchmark blade for the NREL-5MW Reference Wind Turbine, an alternative fixed-speed stall-controlled blade design proposed by Delft University, and a novel alternative design that we proposed for a variable-speed stall-controlled turbine, which offers the potential for more uniform power control and improved annual energy production.. To optimize the power output of the stall-controlled blade we modify the existing designs and study their behavior using the aforementioned aero elastic model.

DYNAMICS AND CONTROL OF REACTIVE DISTILLATION PROCESS FOR MONOMER SYNTHESIS OF POLYCARBONATE PLANTS

Polycarbonate (PC) is an important engineering thermoplastic that is currently produced in large industrial scale using bisphenol A and monomers such as phosgene. Since phosgene is highly toxic, a non-phosgene approach using diphenyl carbonate (DPC) as an alternative monomer, as developed by Asahi Corporation of Japan, is a significantly more environmentally friendly alternative. Other advantages include the use of CO2 instead of CO as raw material and the elimination of major waste water production. However, for the production of DPC to be economically viable, reactive-distillation units are needed to obtain the necessary yields by shifting the reaction-equilibrium to the desired products and separating the products at the point where the equilibrium reaction occurs. In the field of chemical reaction engineering, there are many reactions that are suffering from the low equilibrium constant. The main goal of this research is to determine the optimal process needed to shift the reactions by using appropriate control strategies of the reactive distillation system. An extensive dynamic mathematical model has been developed to help us investigate different control and processing strategies of the reactive distillation units to increase the production of DPC. The high-fidelity dynamic models include extensive thermodynamic and reaction-kinetics models while incorporating the necessary mass and energy balance of the various stages of the reactive distillation units. The study presented in this document shows the possibility of producing DPC via one reactive distillation instead of the conventional two-column, with a production rate of 16.75 tons/h corresponding to start reactants materials of 74.69 tons/h of Phenol and 35.75 tons/h of Dimethyl Carbonate. This represents a threefold increase over the projected production rate given in the literature based on a two-column configuration. In addition, the purity of the DPC produced could reach levels as high as 99.5% with the effective use of controls. These studies are based on simulation done using high-fidelity dynamic models.

Transfer Functions in Artificial Neural Networks - A Simulation-Based Tutorial

Artificial neural networks are based on computational units that resemble basic information processing properties of biological neurons in an abstract and simplified manner. Generally, these formal neurons model an input-output behaviour as it is also often used to characterize biological neurons. The neuron is treated as a black box; spatial extension and temporal dynamics present in biological neurons are most often neglected. Even though artificial neurons are simplified, they can show a variety of input-output relations, depending on the transfer functions they apply. This unit on transfer functions provides an overview of different transfer functions and offers a simulation that visualizes the input-output behaviour of an artificial neuron depending on the specific combination of transfer functions.

Forcing of stratospheric chemistry and dynamics during the Dalton Minimum

The response of atmospheric chemistry and dynamics to volcanic eruptions and to a decrease in solar activity during the Dalton Minimum is investigated with the fully coupled atmosphere–ocean chemistry general circulation model SOCOL-MPIOM (modeling tools for studies of SOlar Climate Ozone Links-Max Planck Institute Ocean Model) covering the time period 1780 to 1840 AD. We carried out several sensitivity ensemble experiments to separate the effects of (i) reduced solar ultra-violet (UV) irradiance, (ii) reduced solar visible and near infrared irradiance, (iii) enhanced galactic cosmic ray intensity as well as less intensive solar energetic proton events and auroral electron precipitation, and (iv) volcanic aerosols. The introduced changes of UV irradiance and volcanic aerosols significantly influence stratospheric dynamics in the early 19th century, whereas changes in the visible part of the spectrum and energetic particles have smaller effects. A reduction of UV irradiance by 15%, which represents the presently discussed highest estimate of UV irradiance change caused by solar activity changes, causes global ozone decrease below the stratopause reaching as much as 8% in the midlatitudes at 5 hPa and a significant stratospheric cooling of up to 2 °C in the mid-stratosphere and to 6 °C in the lower mesosphere. Changes in energetic particle precipitation lead only to minor changes in the yearly averaged temperature fields in the stratosphere. Volcanic aerosols heat the tropical lower stratosphere, allowing more water vapour to enter the tropical stratosphere, which, via HOx reactions, decreases upper stratospheric and mesospheric ozone by roughly 4%. Conversely, heterogeneous chemistry on aerosols reduces stratospheric NOx, leading to a 12% ozone increase in the tropics, whereas a decrease in ozone of up to 5% is found over Antarctica in boreal winter. The linear superposition of the different contributions is not equivalent to the response obtained in a simulation when all forcing factors are applied during the Dalton Minimum (DM) – this effect is especially well visible for NOx/NOy. Thus, this study also shows the non-linear behaviour of the coupled chemistry-climate system. Finally, we conclude that especially UV and volcanic eruptions dominate the changes in the ozone, temperature and dynamics while the NOx field is dominated by the energetic particle precipitation. Visible radiation changes have only very minor effects on both stratospheric dynamics and chemistry.

Reanalysis-driven climate simulation over CORDEX North America domain using the Canadian Regional Climate Model, version 5: model performance evaluation

The performance of reanalysis-driven Canadian Regional Climate Model, version 5 (CRCM5) in reproducing the present climate over the North American COordinated Regional climate Downscaling EXperiment domain for the 1989–2008 period has been assessed in comparison with several observation-based datasets. The model reproduces satisfactorily the near-surface temperature and precipitation characteristics over most part of North America. Coastal and mountainous zones remain problematic: a cold bias (2–6 °C) prevails over Rocky Mountains in summertime and all year-round over Mexico; winter precipitation in mountainous coastal regions is overestimated. The precipitation patterns related to the North American Monsoon are well reproduced, except on its northern limit. The spatial and temporal structure of the Great Plains Low-Level Jet is well reproduced by the model; however, the night-time precipitation maximum in the jet area is underestimated. The performance of CRCM5 was assessed against earlier CRCM versions and other RCMs. CRCM5 is shown to have been substantially improved compared to CRCM3 and CRCM4 in terms of seasonal mean statistics, and to be comparable to other modern RCMs.

Atomic Quantum Simulation of U(N) and SU(N) Abelian Lattice Gauge Theories

Using ultracold alkaline-earth atoms in optical lattices, we construct a quantum simulator for U(N) and SU(N) lattice gauge theories with fermionic matter based on quantum link models. These systems share qualitative features with QCD, including chiral symmetry breaking and restoration at nonzero temperature or baryon density. Unlike classical simulations, a quantum simulator does not suffer from sign problems and can address the corresponding chiral dynamics in real time.

Computer simulation of glioma growth and morphology.

Despite major advances in the study of glioma, the quantitative links between intra-tumor molecular/cellular properties, clinically observable properties such as morphology, and critical tumor behaviors such as growth and invasiveness remain unclear, hampering more effective coupling of tumor physical characteristics with implications for prognosis and therapy. Although molecular biology, histopathology, and radiological imaging are employed in this endeavor, studies are severely challenged by the multitude of different physical scales involved in tumor growth, i.e., from molecular nanoscale to cell microscale and finally to tissue centimeter scale. Consequently, it is often difficult to determine the underlying dynamics across dimensions. New techniques are needed to tackle these issues. Here, we address this multi-scalar problem by employing a novel predictive three-dimensional mathematical and computational model based on first-principle equations (conservation laws of physics) that describe mathematically the diffusion of cell substrates and other processes determining tumor mass growth and invasion. The model uses conserved variables to represent known determinants of glioma behavior, e.g., cell density and oxygen concentration, as well as biological functional relationships and parameters linking phenomena at different scales whose specific forms and values are hypothesized and calculated based on in vitro and in vivo experiments and from histopathology of tissue specimens from human gliomas. This model enables correlation of glioma morphology to tumor growth by quantifying interdependence of tumor mass on the microenvironment (e.g., hypoxia, tissue disruption) and on the cellular phenotypes (e.g., mitosis and apoptosis rates, cell adhesion strength). Once functional relationships between variables and associated parameter values have been informed, e.g., from histopathology or intra-operative analysis, this model can be used for disease diagnosis/prognosis, hypothesis testing, and to guide surgery and therapy. In particular, this tool identifies and quantifies the effects of vascularization and other cell-scale glioma morphological characteristics as predictors of tumor-scale growth and invasion.

Dynamics of a minimal model of interlocked positive and negative feedback loops of transcriptional regulation by cAMP-response element binding proteins.

cAMP-response element binding (CREB) proteins are involved in transcriptional regulation in a number of cellular processes (e.g., neural plasticity and circadian rhythms). The CREB family contains activators and repressors that may interact through positive and negative feedback loops. These loops can be generated by auto- and cross-regulation of expression of CREB proteins, via CRE elements in or near their genes. Experiments suggest that such feedback loops may operate in several systems (e.g., Aplysia and rat). To understand the functional implications of such feedback loops, which are interlocked via cross-regulation of transcription, a minimal model with a positive and negative loop was developed and investigated using bifurcation analysis. Bifurcation analysis revealed diverse nonlinear dynamics (e.g., bistability and oscillations). The stability of steady states or oscillations could be changed by time delays in the synthesis of the activator (CREB1) or the repressor (CREB2). Investigation of stochastic fluctuations due to small numbers of molecules of CREB1 and CREB2 revealed a bimodal distribution of CREB molecules in the bistability region. The robustness of the stable HIGH and LOW states of CREB expression to stochastic noise differs, and a critical number of molecules was required to sustain the HIGH state for days or longer. Increasing positive feedback or decreasing negative feedback also increased the lifetime of the HIGH state, and persistence of this state may correlate with long-term memory formation. A critical number of molecules was also required to sustain robust oscillations of CREB expression. If a steady state was near a deterministic Hopf bifurcation point, stochastic resonance could induce oscillations. This comparative analysis of deterministic and stochastic dynamics not only provides insights into the possible dynamics of CREB regulatory motifs, but also demonstrates a framework for understanding other regulatory processes with similar network architecture.