Stochastic linear scaling for metals and non metals

Total energy electronic structure calculations, based on density functional theory or on the more empirical tight binding approach, are generally believed to scale as the cube of the number of electrons. By using the localisaton property of the high temperature density matrix we present exact deterministic algorithms that scale linearly in one dimension and quadratically in all others. We also introduce a stochastic algorithm that scales linearly with system size. These results hold for metallic and non metallic systems and are substantiated by numerical calculations on model systems.Comment: 9 pages, 2 figure

Evolution of the structure of amorphous ice - from low-density amorphous (LDA) through high-density amorphous (HDA) to very high-density amorphous (VHDA) ice

We report results of molecular dynamics simulations of amorphous ice for pressures up to 22.5 kbar. The high-density amorphous ice (HDA) as prepared by pressure-induced amorphization of Ih ice at T=80 K is annealed to T=170 K at various pressures to allow for relaxation. Upon increase of pressure, relaxed amorphous ice undergoes a pronounced change of structure, ranging from the low-density amorphous ice (LDA) at p=0, through a continuum of HDA states to the limiting very high-density amorphous ice (VHDA) regime above 10 kbar. The main part of the overall structural change takes place within the HDA megabasin, which includes a variety of structures with quite different local and medium-range order as well as network topology and spans a broad range of densities. The VHDA represents the limit to densification by adapting the hydrogen-bonded network topology, without creating interpenetrating networks. The connection between structure and metastability of various forms upon decompression and heating is studied and discussed. We also discuss the analogy with amorphous and crystalline silica. Finally, some conclusions concerning the relation between amorphous ice and supercooled water are drawn.Comment: 11 pages, 12 postscript figures. To be published in The Journal of Chemical Physic

An Efficient and Accurate Car-Parrinello-like Approach to Born-Oppenheimer Molecular Dynamics

We present a new method which combines Car-Parrinello and Born-Oppenheimer molecular dynamics in order to accelerate density functional theory based ab-initio simulations. Depending on the system a gain in efficiency of one to two orders of magnitude has been observed, which allows ab-initio molecular dynamics of much larger time and length scales than previously thought feasible. It will be demonstrated that the dynamics is correctly reproduced and that high accuracy can be maintained throughout for systems ranging from insulators to semiconductors and even to metals in condensed phases. This development considerably extends the scope of ab-initio simulations.Comment: 4 pages, 3 figures; Accepted by Phys. Rev. Lett. for publicatio

Linear scaling electronic structure calculations and accurate sampling with noisy forces

Numerical simulations based on electronic structure calculations are finding ever growing applications in many areas of physics. A major limiting factor is however the cubic scaling of the algorithms used. Building on previous work [F. R. Krajewski and M. Parrinello, Phys.Rev. B71, 233105 (2005)] we introduce a novel statistical method for evaluating the inter-atomic forces which scales linearly with system size and is applicable also to metals. The method is based on exact decomposition of the fermionic determinant and on a mapping onto a field theoretical expression. We solve exactly the problem of sampling the Boltzmann distribution with noisy forces. This novel approach can be used in such diverse fields as quantum chromodynamics, quantum Monte Carlo or colloidal physics.Comment: 5 pages, 2 figure

Polyamorphism of ice at low temperatures from constant-pressure simulations

We report results of MD simulations of amorphous ice in the pressure range 0 - 22.5 kbar. The high-density amorphous ice (HDA) prepared by compression of Ih ice at T = 80 K is annealed to T = 170 K at intermediate pressures in order to generate relaxed states. We confirm the existence of recently observed phenomena, the very high-density amorphous ice and a continuum of HDA forms. We suggest that both phenomena have their origin in the evolution of the network topology of the annealed HDA phase with decreasing volume, resulting at low temperatures in the metastability of a range of densities.Comment: 11 pages, 5 postscript figures. To be published in Physical Review Letter

Predicting crystal structures: the Parrinello-Rahman method revisited

By suitably adapting a recent approach [A. Laio and M. Parrinello, PNAS, 99, 12562 (2002)] we develop a powerful molecular dynamics method for the study of pressure-induced structural transformations. We use the edges of the simulation cell as collective variables. In the space of these variables we define a metadynamics that drives the system away from the local minimum towards a new crystal structure. In contrast to the Parrinello-Rahman method our approach shows no hysteresis and crystal structure transformations can occur at the equilibrium pressure. We illustrate the power of the method by studying the pressure-induced diamond to simple hexagonal phase transition in a model of silicon.Comment: 5 pages, 2 Postscript figures, submitte

Effective reference and current integration for large displacement interface

The most common interface formulations proposed in literature are generally based on the restrictive hypothesis of small strains and small displacements and, even though their application to geometrically nonlinear problems is of paramount interest, only few contributions are available in literature. Motivations are probably due to the difficulties encountered on such formulation, as already mentioned by several authors. A pioneering formulation is the finite displacement three-dimensional interface developed by Ortiz and Pandolfi in [1], where normal and tangential traction components are evaluated with respect to the middle surface in the current configuration, producing a non-symmetric geometric stiffness matrix. More recently, an interface element formulation for geometrical non-linearity and material nonlinearity, which is developed in the reference configuration, has been proposed by Reinoso and Paggi in [2]. The constitutive model is formulated on the local reference, defined by normal axis and tangential axis with respect to the middle surface in the current configuration. The interface formulation generates a non symmetric geometric stiffness matrix, which is simplified by neglecting the non symmetric contribution, in order to reduces computational cost by the use of symmetric solver. The state of the art of cohesive models for the material separation is presented by Mosler and Scheider in [3], focusing the attention on the thermodynamics and variational consistency. In [3] the authors state that many proposed models do not verify fundamental requirements such as thermodynamic principles, frame invariance, equilibrium conditions. Such problems are magnified for anisotropic models in geometrically nonlinear context. Attention is also focused on the unphysical dissipation produced in elastic paths due to unsymmetrical stiffness matrix. Some existing cohesive-zone models are analyzed under conditions of large displacement and large strain by Ottosen et al in [4], and CZMs are also evaluated with respect to thermodynamic consistency and the fundamental laws such as balance of angular momentum and frame invariance. It is shown that in elastic regime only isotropic models, with traction vector aligned to separation displacement vector, fulfill the physical principles, as already shown in [5]. In [6] some cohesive-zone models are compared at finite strain condition, by a wedge test and by a peel test. The paper [6] shows that some models available in literature, or implemented in commercial finite element codes, which integrate the weak form equilibrium condition over the current configuration, produce significant error in terms of fracture energy. On the contrary, models integrated over the reference configuration produce negligible numerical error. The present paper investigates reasons of the different results between current and reference integration schemes. It is shown that interface formulations integrated over current configuration violate energy conservation principle, due to the elastic energy generated by the finite interface elongation with constant elastic stiffness parameters. Moreover, an original mechanical interpretation of the elastic stiffness parameters, defined as a density of elastic springs between the two interface edges, can be considered an effective solution for interface integrated over the current configuration. In fact, the interface elongation modify the density of springs, as well as volume change modifies the mass density, and integration over current configuration and integration over the reference one produce two identical solutions. In the present paper the interface formulation is rigorously developed under large displacement conditions, assuming as local reference for the constitutive model, normal and tangential axes to the middle surface, as already proposed in [1]. The geometric operators in the current configuration, such as the normal and tangential axes to the middle surface and elongation of the middle surface, are defined as functions of nodal displacements, and first order and second order derivatives, with respect to nodal displacements, are developed. Finally, nodal force vector and consistent stiffness matrix are developed for a two-dimensional interface element, showing the symmetry condition of the geometric stiffness matrix, if the second order derivative are not neglected. The proposed interface formulation is implemented in the FEAP finite element code [7] and the cohesive formulation proposed in [8] is considered as constitutive model. Results of numerical some simulations are proposed with times of convergence obtained with a symmetric solver

Accelerating the convergence of path integral dynamics with a generalized Langevin equation

The quantum nature of nuclei plays an important role in the accurate modelling of light atoms such as hydrogen, but it is often neglected in simulations due to the high computational overhead involved. It has recently been shown that zero-point energy effects can be included comparatively cheaply in simulations of harmonic and quasi-harmonic systems by augmenting classical molecular dynamics with a generalized Langevin equation (GLE). Here we describe how a similar approach can be used to accelerate the convergence of path integral (PI) molecular dynamics to the exact quantum mechanical result in more strongly anharmonic systems exhibiting both zero point energy and tunnelling effects. The resulting PI-GLE method is illustrated with applications to a double-well tunnelling problem and to liquid water