Search

Machine learning potentials (MLP) have revolutionized the field of atomistic simulations by describing the atomic interactions with the accuracy of electronic structure methods at a small fraction of the costs. Most current MLPs construct the energy of a system as a sum of atomic energies, which depend on information about the atomic environments provided in form of predefined or learnable feature vectors. If, in addition, non-local phenomena like long-range charge transfer are important, fourth-generation MLPs need to be used, which include a charge equilibration (Qeq) step to take the global structure of the system into account. This Qeq can significantly increase the computational cost and thus can become the computational bottleneck for large systems. In this paper we present a highly efficient formulation of Qeq that does not require the explicit computation of the Coulomb matrix elements resulting in a quasi-linearly scaling method. Moreover, our approach also allows for the efficient calculation of energy derivatives, which explicitly consider the global structure-dependence of the atomic charges as obtained from Qeq. Due to its generality, the method is not restricted to MLPs but can also be applied within a variety of other force fields.

In the dynamic field of materials science, the quest to find optimal structures with low potential energy is of great significance. Over the past two decades, the minima hopping algorithm has emerged as a successful tool in this pursuit. We present a robust, user friendly and efficient implementation of the minima hopping algorithm as a Python library, enhancing in this way the global structure optimization simulations significantly. Our implementation significantly accelerates the exploration the potential energy surfaces, leveraging an MPI parallelization scheme that allows for multi level parallelization. In this scheme, multiple minima hopping processes are running simultaneously communicating their findings to a single database and, therefore, sharing information with each other about which parts of the potential energy surface have already been explored. Also multiple features from several existing implementations such as variable cell shape molecular dynamics and combined atomic position and cell geometry optimization for bulk systems, enhanced temperature feedback and fragmentation fixing for clusters are included in this implementation. Finally, this implementation takes advantage of the Atomic Simulation Environment (ASE) Python library allowing for high flexibility regarding the underlying energy and force evaluation.

This paper presents the results of an extensive structural search of ternary solids containing lutetium, nitrogen and hydrogen. Based on thousands of thermodynamically stable structures, available online, the convex hull of the formation enthalpies is constructed. To obtain the correct energetic ordering, the highly accurate RSCAN DFT functional is used in high quality all-electron calculations. In this way possible pseudopotential errors are eliminated. A novel lutetium hydride structure (HLu$_2$) that is on the convex hull is found in our search. An electron phonon analysis however shows that it is not a candidate structure for near ambient superconductivity. Besides this structure, which appears to have been missed in previous searches, possibly due to different DFT methodologies, our results agree closely with the results of previously published structure search efforts. This shows, that the field of crystal structure prediction has matured to a state where independent methodologies produce consistent and reproducible results, underlining the trustworthiness of modern crystal structure predictions. Hence it is quite unlikely that a structure, that would give rise within standard BCS theory to the superconducting properties, claimed to have been observed by Dasenbrock-Gammon et al. 10.1038/s41586-023-05742-0 , exists. This solidifies the evidence that no structure with conventional superconducting properties exists that could explain the experimental observation made by Dasenbrock-Gammon et al. 10.1038/s41586-023-05742-0

Ground state structures found in nature are in many cases of high symmetry. But structure prediction methods typically render only a small fraction of high symmetry structures. Especially for large crystalline unit cells there are many low energy defect structures. For this reason methods have been developed where either preferentially high symmetry structures are used as input or where the whole structural search is done within a certain symmetry group. In both cases it is necessary to specify the correct symmetry group beforehand. However it can in general not be predicted which symmetry group is the correct one leading to the ground state. For this reason we introduce a potential energy biasing scheme that favors symmetry and where it is not necessary to specify any symmetry group beforehand. On this biased potential energy surface, high symmetry structures will be found much faster than on an unbiased surface and independently of the symmetry group to which they belong. For our two test cases, a $C_{60}$ fullerene and bulk silicon carbide, we get a speedups of 25 and 63. In our data we also find a clear correlation between the similarity of the atomic environments and the energy. In low energy structures all the atoms of a species tend to have similar environments., Comment: 8 pages, 4 figures, 1 table

A fast and reliable geometry optimization algorithm is presented that optimizes atomic positions and lattice vectors simultaneously. Using a series of benchmarks, it is shown that the method presented in this paper outperforms in most cases the standard optimization methods implemented in popular codes such as QUANTUM ESPRESSO and VASP. To motivate the variable cell shape optimization method presented in here, the eigenvalues of the lattice Hessian matrix are investigated thoroughly. It is shown that they change depending on the shape of the cell and the number of particles inside the cell. For certain cell shapes the resulting condition number of the lattice matrix can grow quadratically with respect to the number of particles. By a coordinate transformation which can be applied to all variable cell shape optimization methods, the undesirable conditioning of the lattice Hessian matrix is eliminated.

Atomic fingerprints are commonly used for the characterization of local environments of atoms in machine learning and other contexts. In this work, we study the behavior of two widely used fingerprints, namely the smooth overlap of atomic positions (SOAP) and the atom-centered symmetry functions (ACSF), under finite changes of atomic positions and demonstrate the existence of manifolds of quasi-constant fingerprints. These manifolds are found numerically by following eigenvectors of the sensitivity matrix with quasi-zero eigenvalues. The existence of such manifolds in ACSF and SOAP causes a failure to machine learn four-body interactions such as torsional energies that are part of standard force fields. No such manifolds can be found for the Overlap Matrix (OM) fingerprint due to its intrinsic many-body character., Comment: 8 pages, 5 figures

Machine learning potentials have become an important tool for atomistic simulations in many fields, from chemistry via molecular biology to materials science. Most of the established methods, however, rely on local properties and are thus unable to take global changes in the electronic structure into account, which result from long-range charge transfer or different charge states. In this work we overcome this limitation by introducing a fourth-generation high-dimensional neural network potential that combines a charge equilibration scheme employing environment-dependent atomic electronegativities with accurate atomic energies. The method, which is able to correctly describe global charge distributions in arbitrary systems, yields much improved energies and substantially extends the applicability of modern machine learning potentials. This is demonstrated for a series of systems representing typical scenarios in chemistry and materials science that are incorrectly described by current methods, while the fourth-generation neural network potential is in excellent agreement with electronic structure calculations., Comment: 13 pages, 11 figures

Fingerprint distances, which measure the similarity of atomic environments, are commonly calculated from atomic environment fingerprint vectors. In this work we present the simplex method which can perform the inverse operation, i.e. calculating fingerprint vectors from fingerprint distances. The fingerprint vectors found in this way point to the corners of a simplex. For a large data set of fingerprints, we can find a particular largest volume simplex, whose dimension gives the effective dimension of the fingerprint vector space. We show that the corners of this simplex correspond to landmark environments that can by used in a fully automatic way to analyse structures. In this way we can for instance detect atoms in grain boundaries or on edges of carbon flakes without any human input about the expected environment. By projecting fingerprints on the largest volume simplex we can also obtain fingerprint vectors that are considerably shorter than the original ones but whose information content is not significantly reduced., Comment: 10 pages, 7 figures

Using methods borrowed from machine learning we detect in a fully algorithmic way long range effects on local physical properties in a simple covalent system of carbon atoms. The fact that these long range effects exist for many configurations implies that atomistic simulation methods, such as force fields or modern machine learning schemes, that are based on locality assumptions, are limited in accuracy. We show that the basic driving mechanism for the long range effects is charge transfer. If the charge transfer is known, locality can be recovered for certain quantities such as the band structure energy., Comment: 3 figures

Atomic environment fingerprints are widely used in computational materials science, from machine learning potentials to the quantification of similarities between atomic configurations. Many approaches to the construction of such fingerprints, also called structural descriptors, have been proposed. In this work, we compare the performance of fingerprints based on the Overlap Matrix(OM), the Smooth Overlap of Atomic Positions (SOAP), Behler-Parrinello atom-centered symmetry functions (ACSF), modified Behler-Parrinello symmetry functions (MBSF) used in the ANI-1ccx potential and the Faber-Christensen-Huang-Lilienfeld (FCHL) fingerprint under various aspects. We study their ability to resolve differences in local environments and in particular examine whether there are certain atomic movements that leave the fingerprints exactly or nearly invariant. For this purpose, we introduce a sensitivity matrix whose eigenvalues quantify the effect of atomic displacement modes on the fingerprint. Further, we check whether these displacements correlate with the variation of localized physical quantities such as forces. Finally, we extend our examination to the correlation between molecular fingerprints obtained from the atomic fingerprints and global quantities of entire molecules., Comment: 18 pages, 11 figures

Monte Carlo simulations are a powerful tool to investigate the thermodynamic properties of atomic systems. In practice however, sampling of the complete configuration space is often hindered by high energy barriers between different regions of configuration space which can make ergodic sampling completely infeasible within accessible simulation times. Although several extensions to the conventional Monte Carlo scheme have been developed, that enable the treatment of such systems, these extensions often entail substantial computational cost or rely on the harmonic approximation. In this work we propose an exact method called Funnel Hopping Monte Carlo (FHMC) that is inspired by the the ideas of smart darting but is more efficient. Gaussian mixtures are used to approximate the Boltzmann distribution around local energy minima which are then used to propose high quality Monte Carlo moves that enable the Monte Carlo simulation to directly jump between different funnels. To fit the Gaussian mixtures we developed an extended version of the expectation-maximization algorithm that is able to take advantage of the high symmetry present in many low energy configurations. We demonstrate the methods performance on the example of the 38 as well as the 75 atom Lennard-Jones clusters which are well known for their double funnel energy landscapes that prevent ergodic sampling with conventional Monte Carlo simulations. By integrating FHMC into the parallel tempering scheme we were able to reduce the number of steps required until convergence of the simulation significantly., Comment: 12 pages, 6 figures

Finding complex reaction and transformation pathways, involving many intermediate states, is in general not possible on the DFT level with existing simulation methods due to the very large number of required energy and force evaluations. This is due to a large extent to the fact that for complex reactions, it is not possible to determine which atom in the educt is mapped onto which atom in the product. Trying out all possible atomic index mappings is not feasible because of the factorial increase in the number of possible mappings. By using a penalty function that is invariant under index permutations, we can bias the potential energy surface in such a way that it obtains the characteristics of a structure seeker whose global minimum is the product. By performing a Minima Hopping based global optimization on this biased potential energy surface we can rapidly find intermediate states that lead into the global minimum. Based on this information we can then extract the full reaction pathway. We first demonstrate for a benchmark system, namely LJ38 that our method allows to preferentially find intermediate states that are relevant for the lowest energy reaction pathway and that we, therefore, need a much smaller number of intermediate states than previous methods to find the lowest energy reaction pathway. Finally, we apply the method to two real systems, C60 and C20H20 and show that the found reaction pathway contains valuable information on how the system can be synthesized.

Using multiwavelets, we have obtained total energies and corresponding atomization energies for the GGA-PBE and hybrid-PBE0 density functionals for a test set of 211 molecules with an unprecedented and guaranteed $\mu$Hartree accuracy. These quasi-exact references allow us to quantify the accuracy of standard all-electron basis sets that are believed to be highly accurate for molecules, such as Gaussian-type orbitals (GTOs), all-electron numeric atom-centered orbitals (NAOs) and full-potential augmented plane wave (APW) methods. We show that NAOs are able to achieve the so-called chemical accuracy (1 kcal/mol) for the typical basis set sizes used in applications, for both total and atomization energies. For GTOs, a triple-zeta quality basis has mean errors of ~10kcal/mol in total energies, while chemical accuracy is almost reached for a quintuple-zeta basis. Due to systematic error cancellations, atomization energy errors are reduced by almost an order of magnitude, placing chemical accuracy within reach also for medium to large GTO bases, albeit with significant outliers. In order to check the accuracy of the computed densities, we have also investigated the dipole moments, where in general, only the largest NAO and GTO bases are able to yield errors below 0.01 Debye. The observed errors are similar across the different functionals considered here.

Measuring similarities/dissimilarities between atomic structures is important for the exploration of potential energy landscapes. However, the cell vectors together with the coordinates of the atoms, which are generally used to describe periodic systems, are quantities not suitable as fingerprints to distinguish structures. Based on a characterization of the local environment of all atoms in a cell we introduce crystal fingerprints that can be calculated easily and allow to define configurational distances between crystalline structures that satisfy the mathematical properties of a metric. This distance between two configurations is a measure of their similarity/dissimilarity and it allows in particular to distinguish structures. The new method is an useful tool within various energy landscape exploration schemes, such as minima hopping, random search, swarm intelligence algorithms and high-throughput screenings.

Optimizations of atomic positions belong to the most commonly performed tasks in electronic structure calculations. Many simulations like global minimum searches or characterizations of chemical reactions require performing hundreds or thousands of minimizations or saddle computations. To automatize these tasks, optimization algorithms must not only be efficient, but also very reliable. Unfortunately computational noise in forces and energies is inherent to electronic structure codes. This computational noise poses a sever problem to the stability of efficient optimization methods like the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm. We here present a technique that allows obtaining significant curvature information of noisy potential energy surfaces. We use this technique to construct both, a stabilized quasi-Newton minimization method and a stabilized quasi-Newton saddle finding approach. We demonstrate with the help of benchmarks that both the minimizer and the saddle finding approach are superior to comparable existing methods.

The Minima Hopping global optimization method uses physically realizable molecular dynamics moves in combination with an energy feedback that guarantees the escape from any potential energy funnel. For the purpose of finding reactions pathways, we argue that Minima Hopping is particularly suitable as a guide through the potential energy landscape and as a generator for pairs of minima that can be used as input structures for methods capable of finding transition states between two minima. For Lennard-Jones benchmark systems we compared this Minima Hopping guided path search method to a known approach for the exploration of potential energy landscapes that is based on deterministic mode-following. Although we used a stabilized mode-following technique that reliably allows to follow distinct directions when escaping from a local minimum, we observed that Minima Hopping guided path search is far superior in finding lowest-barrier reaction pathways. We therefore suggest that Minima Hopping guided path search can be used as a simple and efficient way to identify energetically low-lying chemical reaction pathways. Finally we applied the Minima Hopping guided path search approach to 75-atom and 102-atom Lennard Jones systems. For the 75-atom system we found pathways whose highest energies are significantly lower than the highest energy along the previously published lowest-barrier pathway. Furthermore, many of these pathways contain a smaller number of intermediate transition states than the previously publish lowest-barrier pathway. In case of the 102-atom system Minima Hopping guided path search found a previously unknown and energetically low-lying funnel., Comment: 15 pages, 8 figures

A 3-dimensional GPU Poisson solver is developed for all possible combinations of free and periodic boundary conditions (BCs) along the three directions. It is benchmarked for various grid sizes and different BCs and a significant performance gain is observed for problems including one or more free BCs. The GPU Poisson solver is also benchmarked against two different CPU implementations of the same method and a significant amount of acceleration of the computation is observed with the GPU version., Comment: 10 pages, 5 figures

A novel allotrope of carbon with $P2/m$ symmetry was identified during an \emph{ab-initio} minima-hopping structural search which we call $M10$-carbon. This structure is predicted to be more stable than graphite at pressures above 14.4 GPa and consists purely of $sp^3$ bonds. It has a high bulk modulus and is almost as hard as diamond. A comparison of the simulated X-ray diffraction pattern shows a good agreement with experimental results from cold compressed graphite., Comment: 3 pages, 3 figures

Through a systematic structural search we found an allotrope of carbon with Cmmm symmetry which we predict to be more stable than graphite for pressures above 10 GPa. This material, which we refer to as Z-carbon, is formed by pure sp3 bonds and is the only carbon allotrope which provides an excellent match to unexplained features in experimental X-ray diffraction and Raman spectra of graphite under pressure. The transition from graphite to Z-carbon can occur through simple sliding and buckling of graphene sheets. Our calculations predict that Z-carbon is a transparent wide band gap semiconductor with a hardness comparable to diamond., Comment: 4 pages, 5 figures

A structure prediction method is presented based on the Minima Hopping method. Optimized moves on the configurational enthalpy surface are performed to escape local minima using variable cell shape molecular dynamics by aligning the initial atomic and cell velocities to low curvature directions of the current minimum. The method is applied to both silicon crystals and binary Lennard-Jones mixtures and the results are compared to previous investigations. It is shown that a high success rate is achieved and a reliable prediction of unknown ground state structures is possible., Comment: 9 pages, 6 figures, novel approach in structure prediction, submitted to the Journal of Chemical Physics

Endohedrally doped Si20 fullerenes appear as appealing building blocks for nanoscale materials. We investigate their structural stability with an unbiased and systematic global geometry optimization method within density-functional theory. For a wide range of metal doping atoms, it was sufficient to explore the Born Oppenheimer surface for only a moderate number of local minima to find structures that clearly differ from the initial endohedral cages, but are considerably more favorable in terms of energy. Previously proposed structures are thus all metastable., Comment: 4 pages, 1 figure

The accuracy of the energy landscape of silicon systems obtained from various density functional methods, a tight binding scheme and force fields is studied. Quantum Monte Carlo results serve as quasi exact reference values. In addition to the well known accuracy of DFT methods for geometric ground states and metastable configurations we find that DFT methods give a similar accuracy for transition states and thus a good overall description of the energy landscape. On the other hand, force fields give a very poor description of the landscape that are in most cases too rugged and contain many fake local minima and saddle points or ones that have the wrong height.

The implementation of a full electronic structure calculation code on a hybrid parallel architecture with Graphic Processing Units (GPU) is presented. The code which is on the basis of our implementation is a GNU-GPL code based on Daubechies wavelets. It shows very good performances, systematic convergence properties and an excellent efficiency on parallel computers. Our GPU-based acceleration fully preserves all these properties. In particular, the code is able to run on many cores which may or may not have a GPU associated. It is thus able to run on parallel and massive parallel hybrid environment, also with a non-homogeneous ratio CPU/GPU. With double precision calculations, we may achieve considerable speedup, between a factor of 20 for some operations and a factor of 6 for the whole DFT code., Comment: 14 pages, 8 figures

We present a fast and accurate method to calculate the electrostatic energy and forces of interacting particles with the boundary conditions appropriate to surfaces, i.e periodic in the two directions parallel to the surface and free in the perpendicular direction. In the spirit of the Ewald method the problem is divided into a short range and long range part. The charge density responsible for the long range part is represented by plane waves in the periodic directions and by finite elements in the non-periodic direction. Our method has computational complexity of O(N_g log(N_g)) with a very small prefactor, where N_g is the number of grid points., Comment: It has 19 pages and 4 figures

The Bell-Evans-Polanyi principle that is valid for a chemical reaction that proceeds along the reaction coordinate over the transition state is extended to molecular dynamics trajectories that in general do not cross the dividing surface between the initial and the final local minima at the exact transition state. Our molecular dynamics Bell-Evans-Polanyi principle states that low energy molecular dynamics trajectories are more likely to lead into the basin of attraction of a low energy local minimum than high energy trajectories. In the context of global optimization schemes based on molecular dynamics our molecular dynamics Bell-Evans-Polanyi principle implies that using low energy trajectories one needs to visit a smaller number of distinguishable local minima before finding the global minimum than when using high energy trajectories.

An algorithm for fast calculation of the Coulombic forces and energies of point particles with free boundary conditions is proposed. Its calculation time scales as N log N for N particles. This novel method has lower crossover point with the full O(N^2) direct summation than the Fast Multipole Method. The forces obtained by our algorithm are analytical derivatives of the energy which guarantees energy conservation during a molecular dynamics simulation. Our algorithm is very simple. An MPI parallelised version of the code can be downloaded under the GNU General Public License from the website of our group., Comment: 19 pages, 11 figures, submitted to: Journal of Chemical Physics

Wavelets are a powerful new mathematical tool which offers the possibility to treat in a natural way quantities characterized by several length scales. In this article we will show how wavelets can be used to solve partial differential equations which exhibit widely varying length scales and which are therefore hardly accessible by other numerical methods. As a benchmark calculation we solve Poisson's equation for a 3-dimensional Uranium dimer. The length scales of the charge distribution vary by 4 orders of magnitude in this case. Using lifted interpolating wavelets the number of iterations is independent of the maximal resolution and the computational effort therefore scales strictly linearly with respect to the size of the system.