The Royal Society RAMP modelling initiative.

Normally, science proceeds following a well-established set of principles. Studies are done with an emphasis on correctness, are submitted to a journal editor who evaluates their relevance, and then undergo anonymous peer review by experts before publication in a journal and acceptance by the scientific community via the open literature. This process is slow, but its accuracy has served all fields of science well. In an emergency situation, different priorities come to the fore. Research and review need to be conducted quickly, and the target audience consists of policymakers. Scientists must jostle for the attention of non-specialists without sacrificing rigour, and must deal not only with peer assessment but also with media scrutiny by journalists who may have agendas other than ensuring scientific correctness. Here, we describe how the Royal Society coordinated efforts of diverse scientists to help model the coronavirus epidemic. This article is part of the theme issue 'Technical challenges of modelling real-life epidemics and examples of overcoming these'.

Technical challenges of modelling real-life epidemics and examples of overcoming these.

The coronavirus disease 2019 (COVID-19) pandemic has highlighted the importance of mathematical modelling in informing and advising policy decision-making. Effective practice of mathematical modelling has challenges. These can be around the technical modelling framework and how different techniques are combined, the appropriate use of mathematical formalisms or computational languages to accurately capture the intended mechanism or process being studied, in transparency and robustness of models and numerical code, in simulating the appropriate scenarios via explicitly identifying underlying assumptions about the process in nature and simplifying approximations to facilitate modelling, in correctly quantifying the uncertainty of the model parameters and projections, in taking into account the variable quality of data sources, and applying established software engineering practices to avoid duplication of effort and ensure reproducibility of numerical results. Via a collection of 16 technical papers, this special issue aims to address some of these challenges alongside showcasing the usefulness of modelling as applied in this pandemic. This article is part of the theme issue 'Technical challenges of modelling real-life epidemics and examples of overcoming these'.

Phase behavior of the quantum Lennard-Jones solid.

The Lennard-Jones (LJ) potential is perhaps one of the most widely used models for the interaction of uncharged particles, such as noble gas solids. The phase diagram of the classical LJ solid is known to exhibit transitions between hcp and fcc phases. However, the phase behavior of the quantum LJ solid remains unknown. Thermodynamic integration based on path integral molecular dynamics (PIMD) and lattice dynamics calculations are used to study the phase stability of the hcp and fcc LJ solids. The hcp phase is shown to be stabilized by quantum effects in PIMD, while fcc is shown to be favored by lattice dynamics, which suggests a possible re-entrant low pressure fcc phase for highly quantum systems. Implications for the phase stability of noble gas solids are discussed. For parameters equating to helium, the expansion due to zero-point vibrations is associated with quantum melting: neither crystal structure is stable at zero pressure.

Development of an interatomic potential for the simulation of defects, plasticity, and phase transformations in titanium.

New interatomic potentials describing defects, plasticity, and high temperature phase transitions for Ti are presented. Fitting the martensitic hcp-bcc phase transformation temperature requires an efficient and accurate method to determine it. We apply a molecular dynamics method based on determination of the melting temperature of competing solid phases, and Gibbs-Helmholtz integration, and a lattice-switch Monte Carlo method: these agree on the hcp-bcc transformation temperatures to within 2 K. We were able to develop embedded atom potentials which give a good fit to either low or high temperature data, but not both. The first developed potential (Ti1) reproduces the hcp-bcc transformation and melting temperatures and is suitable for the simulation of phase transitions and bcc Ti. Two other potentials (Ti2 and Ti3) correctly describe defect properties and can be used to simulate plasticity or radiation damage in hcp Ti. The fact that a single embedded atom method potential cannot describe both low and high temperature phases may be attributed to neglect of electronic degrees of freedom, notably bcc has a much higher electronic entropy. A temperature-dependent potential obtained from the combination of potentials Ti1 and Ti2 may be used to simulate Ti properties at any temperature.

Paramagnetic and glass transitions in sudoku.

We study the statistical mechanics of a model glassy system based on sudoku, a familiar and popular mathematical puzzle. Sudoku puzzles provide a very rare example of a class of frustrated systems with a unique ground state without symmetry. Here, the puzzle is recast as a thermodynamic system where the number of violated rules defines the energy. We use Monte Carlo simulation to show that the "sudoku Hamiltonian" exhibits two transitions as a function of temperature, a paramagnetic, and a glass transition. Of these, the intermediate condensed phase is the only one that visits the ground state (i.e., it solves the puzzle, though this is not the purpose of the study). Both transitions are associated with an entropy change, paramagnetism measured from the dynamics of the Monte Carlo run, showing a peak in specific heat, while the residual glass entropy is determined by finding multiple instances of the glass by repeated annealing. There are relatively few such simple models for frustrated or glassy systems that exhibit both ordering and glass transitions; sudoku puzzles are unique for the ease with which they can be obtained, with the proof of the existence of a unique ground state via the satisfiability of all constraints. Simulations suggest that in the glass phase there is an increase in information entropy with lowering temperature. In fact, we have shown that sudoku puzzles have the type of rugged energy landscape with multiple minima that typifies glasses in many physical systems. This puzzling result is a manifestation of the paradox of the residual glass entropy. These readily available puzzles can now be used as solvable model Hamiltonian systems for studying the glass transition.

Europium-IV: an incommensurately modulated crystal structure in the lanthanides.

High-resolution x-ray powder-diffraction experiments were performed on europium metal at high pressure up to 50 GPa. At variance with previous reports, the hcp phase of Eu was observed to be stable not only to 18 GPa, but to 31.5 GPa. At 31.5(5) GPa, europium transforms to a phase (Eu-IV) with an incommensurately modulated monoclinic crystal structure with superspace group C2/c(q(1)0q(3))00. This new phase was observed to be stable to ~37.0 GPa, where another phase transition was observed. Eu-IV is the first phase in the lanthanide elements with an incommensurate crystal structure.

Lattice dynamics of dense lithium.

We report low-frequency high-resolution Raman spectroscopy and ab-initio calculations on dense lithium from 40 to 200 GPa at low temperatures. Our experimental results reveal rich first-order Raman activity in the metallic and semiconducting phases of lithium. The computed Raman frequencies are in excellent agreement with the measurements. Free energy calculations provide a quantitative description and physical explanation of the experimental phase diagram only when vibrational effect are correctly treated. The study underlines the importance of zero-point energy in determining the phase stability of compressed lithium.

Crystal structures of dense lithium: a metal-semiconductor-metal transition.

Ab initio random structure searching and single-crystal x-ray diffraction have been used to determine the full structures of three phases of lithium, recently discovered at low temperature above 60 GPa. A structure with C2mb symmetry, calculated to be a poor metal, is proposed for the oC88 phase (60-65 GPa). The oC40 phase (65-95 GPa) is found to have a lowest-enthalpy structure with C2cb symmetry, in excellent agreement with the x-ray data. It is calculated to be a semiconductor with a band gap of ∼1 eV at 90 GPa. oC24, stable above 95 GPa, has the space group Cmca, and refined atomic coordinates are in excellent agreement with previous calculations.

Potassium under pressure: a pseudobinary ionic compound.

Experimentally, we have found that among the "complicated" phases of potassium at intermediate pressures is one which has the same space group as the double hexagonal-close-packed structure, although its atomic coordination is completely different. Calculations on this P6(3)/mmc (hP4) structure as a function of pressure show three isostructural transitions and three distinctive types of chemical bonding: free electron, ionic, and metallic. Interestingly, relationships between localized metallic structures and ionic compounds are found.

Lattice-switch Monte Carlo simulation for binary hard-sphere crystals.

We show how to generalize the lattice-switch Monte Carlo method to calculate the phase diagram of a binary system. A global coordinate transformation is combined with a modification of particle diameters, enabling the multicomponent system in question to be explored and directly compared to a suitable reference state in a single Monte Carlo simulation. We use the method to evaluate the free energies of binary hard-sphere crystals. Calculations at moderate size ratios alpha=0.58 and 0.73 are in agreement with previous results, and confirm AB2 and AB13 as stable structures. We also find that the AB(CsCl) structure is not entropically stable at the size ratio and volume where it has been reported experimentally, and therefore that those observations cannot be explained by packing effects alone.

Constrained dynamics and extraction of normal modes from ab initio molecular dynamics: application to ammonia.

We present a methodology for extracting phonon data from ab initio Born-Oppenheimer molecular dynamics calculations of molecular crystals. Conventional ab initio phonon methods based on perturbations are difficult to apply to lattice modes because the perturbation energy is dominated by intramolecular modes. We use constrained molecular dynamics to eliminate the effect of bond bends and stretches and then show how trajectories can be used to isolate and define in particular, the eigenvalues and eigenvectors of modes irrespective of their symmetry or wave vector. This is done by k-point and frequency filtering and projection onto plane wave states. The method is applied to crystalline ammonia: the constrained molecular dynamics allows a significant speed-up without affecting structural or vibrational modes. All Gamma point lattice modes are isolated: the frequencies are in agreement with previous studies; however, the mode assignments are different.

Magnetically induced immiscibility in the Ising model of FeCr stainless steel.

The iron-chromium alloy system has an unexplained anomaly: although there is a broad miscibility gap it appears to be favorable for chromium to dissolve in iron. This is consistent with ab initio calculation, but no simpler physically intuitive picture has been presented. Here it is shown that the Ising model, based on the bcc lattice with antiferromagnetic and ferromagnetic species, has the potential to exhibit similar behavior, with a skew miscibility gap arising from the solubility of antiferromagnetic species on nonadjacent sites. Essential characteristics of stainless steel (high Cr solubility and surface segregation) are correctly reproduced. Under some conditions, magnetization increases with temperature. The equilibrium miscibility gap due to mixed magnetism and segregation-driven positive dM/dT are fundamental features of the bcc Ising model itself, not just FeCr.

Mutation of albedo and growth response produces oscillations in a spatial Daisyworld.

We extended a two-dimensional cellular automaton (CA) Daisyworld to include mutation of optimal growth temperature as well as mutation of albedo. Thus, the organisms (daisies) can adapt to prevailing environmental conditions or evolve to alter their environment. We find the resulting system oscillates with a period of hundreds of daisy generations. Weaker and less regular oscillations exist in previous daisyworld models, but they become much stronger and more regular here with mutation in the growth response. Despite the existence of a particular combination of mean albedo and optimum individual growth temperature which maximises growth, we find that this global state is unstable with respect to mutations which lower absolute growth rate, but increase marginal growth rate. The resulting system oscillates with a period that is found to decrease with increasing death rate, and to increase with increasing heat diffusion and heat capacity. We speculate that the origin of this oscillation is a Hopf bifurcation, previously predicted in a zero-dimensional system.

Microscopic model of diffusion limited aggregation and electrodeposition in the presence of leveling molecules.

We present a cellular automata approach for microscopic modeling of the effect of unbinding in diffusion limited aggregation. The automata represent active particles, which are able to change their internal state and affect their neighbors. The geometry resembles electrochemical deposition-"ions" diffuse at random from the top of a container until encountering an aggregate in contact with the bottom, to which they stick. This model exhibits dendritic (fractal) growth in the diffusion limited case. The addition of a field eliminates the fractal nature but the density remains low. The addition of molecules that unbind atoms from the aggregate transforms the deposit to a 100% dense one (in three dimensions). The molecules are remarkably adept at avoiding being trapped. This mimics the effect of so-called leveler molecules which are used in electrochemical deposition.

Stabilization of large generalized Lotka-Volterra foodwebs by evolutionary feedback.

Conventional ecological models show that complexity destabilizes foodwebs, suggesting that foodwebs should have neither large numbers of species nor a large number of interactions. However, in nature the opposite appears to be the case. Here we show that if the interactions between species are allowed to evolve within a generalized Lotka-Volterra model such stabilizing feedbacks and weak interactions emerge automatically. Moreover, we show that trophic levels also emerge spontaneously from the evolutionary approach, and the efficiency of the unperturbed ecosystem increases with time. The key to stability in large foodwebs appears to arise not from complexity per se but from evolution at the level of the ecosystem which favors stabilizing (negative) feedbacks.

Maximization principles and daisyworld.

We investigate whether the equilibrium time-averaged state of a self-organizing system with many internal degrees of freedom, 2D-daisyworld, can be described by optimizing a single quantity. Unlike physical systems where a principle of maximum energy production has been observed, daisyworld follows evolutionary dynamics rather than Hamiltonian dynamics. We find that this is sufficient to invalidate the maximum entropy production principle, finding instead a different principle, that the system self-organizes to a state which maximizes the amount of life.

Devil's staircase in kinetically limited growth.

The devil's staircase is a term used to describe surface or an equilibrium phase diagram in which various ordered facets or phases are infinitely closely packed as a function of some model parameter. A classic example is a one-dimensional Ising model [P. Bak and R. Bruinsma, Phys. Rev. Lett. 49, 249 (1982)] wherein long-range and short-range forces compete, and the periodicity of the gaps between minority species covers all rational values. In many physical cases, crystal growth proceeds by adding surface layers that have the lowest energy, but are then frozen in place. The emerging layered structure is not the thermodynamic ground state, but is uniquely defined by the growth kinetics. It is shown that for such a system, the grown structure tends to the equilibrium ground state via a devil's staircase traversing an infinity of intermediate phases. It would be extremely difficult to deduce the simple growth law based on measurement made on such a grown structure.

Lattice-switch Monte Carlo method: application to soft potentials.

The lattice-switch Monte Carlo method, recently introduced and applied in the context of hard spheres, is extended to particles interacting through a soft potential. The method utilizes a transformation that switches between configurations of two different crystalline structures, allowing the phase space of both structures to be explored in a single simulation and the difference between their free energies to be determined directly. We apply the method to determine the fcc-hcp crystalline phase behavior of the classical Lennard-Jones solid.

Pack formation in cycling and orienteering.

In cycling and orienteering competitions, competitors can become bunched into packs, which may mask an individual's true ability. Here we model this process with a view to determining when competitors' times are determined more by others than by their own ability. Our results may prove useful in helping to stage events so that pack formation can be avoided.