link to homepage

Institute for Advanced Simulation (IAS)

Navigation and service

Forschungszentrum Jülich - Workshop HYBRID2013, List of Abstracts

List of Abstracts


Tersoff-Abell type bond-order potentials: A review
Karsten Albe

Bond-order potentials of the Tersoff-Abell type were originally proposed for modelling covalent materials, like carbon and silicon and later extended to other materials including various compounds.

In the contribution I will give a brief introduction into the formalism and discuss fitting strategies as well as various parametrizations for modelling compound semiconductors, magnetic and non-magnetic metals. An extension by a simplified charge-transfer scheme will be presented and several examples for successes and shortcomings of this type of bond-order potential will be given.

Facilitating the development, selection, and use of interatomic potentials (force fields) through the development of robust tools and critical comparisons
Zachary T. Trautt and Chandler A. Becker

In the field of interatomic potentials, a researcher is often in the situation of feast or famine. In the first case, there are many interatomic potentials available for what is nominally the same material, and a user must somehow select one. In the second case, there may be only a few options, none of which is necessarily good enough for the intended application.

We will discuss our efforts to help users select and evaluate interatomic potentials through the development of property calculation methods and their application to the interatomic potentials available through the NIST Interatomic Potentials Repository. These methods are implemented as Python scripts and IPython notebooks in order to facilitate reproducibility and increased scientific openness. Additionally, we will discuss our efforts to compile and use critically assessed reference data for robust comparisons and decision-making.

Sequential multiscale simulations with potfit: First principles data for classical molecular dynamics
Peter Brommer

Download talk as PDF.

Large-scale molecular dynamics (MD) simulations are only feasible with classical interaction potentials or force fields, where all electronic degrees of freedom are summarised in an effective interaction between nuclei. This now only depends on the positions of the atoms and can be evaluated rapidly and efficiently. Unfortunately, determining the parameters of such an effective potential is far from trivial. Traditionally, they were chosen to best reproduce a number of experimentally determined quantities of interest, like elastic constants. However, there is a disconcerting lack of usable information, particularly for more complex materials, where not even the atomic positions are sufficiently well known. Force matching, developed by Ercolessi and Adams about 20 years ago [1], circumvents this by directly using forces on individual atoms determined with first-principles methods as reference data, implying that if a potential gets the forces correct, it will reproduce the atomic dynamics correctly and thus all quantities that are derived from that. The open-source tool potfit is an implementation of the force matching method [2]. Originally developed to determine interaction potentials for quasicrystals and other structurally complex metallic alloys, it has since been extended to cover various potential models for metals and ionic solids.

The program is available from In this talk, I'll discuss the main features of the program as well as recent advances in using potfit for electron-temperature dependent potentials.

[1] Ercolessi and Adams, Europhys. Lett., 1994, 26, 583.
[2] Brommer and Gähler, Modell. Simul. Mater. Sci. Eng., 2007, 15, 295–304.

Force-field based modeling of redox reactions
Wolf Dapp

Download talk as PDF.

Calculating the Coulomb forces in Molecular Dynamics (MD), the effective atomic charge on each atom must be known accurately. While there are several approaches (QE, AACT) for this task, they each have shortcomings. Split Charge Equilibration (SQE) is a method combining the beneficial aspects of both QE and AACT, while avoiding their deficiencies, such as predicting a wrong charge distribution among dissociation of dimers.

To date, partial atomic charges are either assumed constant or assigned according to minimization principles that are smooth and unique functions of the instantaneous atomic coordinates. Thus, charge-transfer history effects cannot be accounted for. Such occur for instance during contact-induced charge transfer between two solid bodies. In particular, redox reactions involve a quasi-discontinuous change of the electronic state, without significant atomic rearrangement.

RedoxSQE, an extension to SQE, assigns discrete oxidation states to each atom. This enables us to model redox reactions. This talk will show proof-of-concept applications, including a complete atomistic nano battery, which self-consistently shows many generic properties of macroscopic batteries. Among other applications, redoxSQE can help to better understand the processes happening at the electrode-electrolyte interface of a real battery.

Bond-Order Potentials: From the electronic structure to million atom simulations
Ralf Drautz

The analytic Bond-Order Potentials (BOPs) are derived from density functional theory (DFT) by two coarse-graining steps. In a first step the electronic structure is simplified to the tight-binding bond approximation by expanding the Hohenberg-Kohn-Sham energy functional to second order in charge and magnetic fluctuations. The second step achieves a local representation of the energy by rewriting the tight-binding approximation in a moments expansion. The BOPs encompass covalent bond formation, charge transfer and non-collinear magnetism within their derived functional form. The analytic BOP framework provides efficient and exact analytic forces and torques that enable simulations with millions of atoms. The BOP parameters follow in analogy by downfolding of the DFT wave functions on a optimized basis and interpolation of the resulting matrix elements.

In this talk I will outline the derivation of the BOPs, discuss their parameterization and give examples of the applications of the BOPs to modelling metals and steels.

Study of grain boundary deformation mechanisms in cemented carbides using a model potential for the W-C-Co system
Martin Gren, Martin Petisme, and Göran Wahnström

Download talk as PDF.

Cemented carbide is a composite material consisting of tungsten carbide grains cemented in a metallic binder. Cobalt is usually the metallic binder. A cemented carbide is produced by means of powder metallurgy in which WC and Co powders are mixed, pressed, and then sintered together. The material inherits significant hardness from the continuous skeleton of hard WC grains while the ductile binder phase serves to increase the toughness. These properties make cemented carbide ubiquitous in a multitude of applications, such as high-speed machining of steels, rock drilling, and for wear parts.

As coatings have increased the material's resistance to abrasive wear, the plastic deformation of the bulk material is often what limits the tool life in practice at high cutting speeds. Grain boundary sliding (GBS) is thought to be the dominating deformation mechanism in WC–Co based cemented carbides at temperatures approaching the Co binder phase melting point.

In order to study the effect of Co infiltration of WC/WC grain boundaries on the sliding, we have performed GBS simulations for bicrystal systems using molecular dynamics. To get good statistical properties and for geometrical reasons the number of atoms needed to model realistic systems lie in the thousands to millions. Also, the number of time steps needed to simulate sliding for relevant sliding distances lie in the millions. This makes first-principles methods such as density functional theory currently impractical. We have therefore used a classical inter-atomic potential for the W-C-Co system. The potential is of the Tersoff form. The interactions for the W-C system are taken from a W-C-H potential found in literature. We have fitted the remaining interactions involving Co using force matching to describe dissolved W and C in Co, and WC/Co interfaces. The Large-scale Atomic/Molecular Massively Parallel Simulator (lammps) software was used to perform the molecular dynamics simulations.

It is experimentally known that a majority of the grain boundaries involve low-energy prismatic or basal planes. For this reason we chose to study sliding of two model grain boundary systems, one with a grain boundary with coincidence index 2 including a basal plane, and one with a grain boundary with coincidence index 4 including a prismatic plane.

We have studied the effect on sliding of segregated Co atoms and thin Co films of varying thickness at 500 K, 1000 K, 1500 K, and 2000 K, where 2000K is above the melting point for Cobalt. The shear stresses for constant sliding rates have been compared and it is concluded that GBS is significantly facilitated by nanometer-thick Co films. Further, we observe that submonolayer segregation of Co generally increases the resistance towards sliding. However, for the grain boundary with coincidence index 2 at temperatures of 1000 K and below this effect is not observed. We found this to be due to the fact that the 0.5 monolayer of Co film in theses cases remains ordered during the sliding process. At higher temperatures for this grain boundary and at all studied temperatures for the other grain boundary the film instead becomes disordered, which hinders the sliding.

Recent developments and simulations utilizing bond-order potentials
Judith A. Harrison

Covalent materials, e.g., Si and C, form strong directional bonds. This poses a challenge for potential development for this important class of materials. Tersoff was the first to attempt to incorporate the structural chemistry of covalently bonded systems into an empirical potential energy function. The development of bond-order potentials (BOP) for C, Si, and Ge has spawned other potentials that utilize the bond-order approach [Tersoff, Phys. Rev. B, 39 (1989) 5566]. Using Tersoff’s BOP for C, Brenner created the reactive empirical bond-order (REBO) potential for hydrocarbons and used it to model diamond growth [Brenner, Phys. Rev. B, 42 (1990) 9458]. The first attempts to construct bond-order potentials capable of modeling 3 different atom types were two independently developed C-Si-H potentials, based on the REBO formalism [Beardmore & Smith, Phil. Mag. A, 74 (1996) 1439; Dyson & Smith, Surf. Sci., 355 (1996) 140]. These potentials inherit the strengths and weaknesses of the underlying potentials. Because the REBO potential is unable to reproduce the elastic constants of diamond and graphite, the C-Si-H potentials also have poor elastic properties. Thus, the functional form of the REBO was updated to create the second-generation REBO, or REBO2, potential [Brenner, et al., J. Phys. C., 14 (2002) 783]. The functional form and expanded fitting database resulted in significantly better descriptions of bond energies, lengths, force constants for hydrocarbon molecules. In addition, the zero-Kelvin elastic properties, interstitial defect energies, and surface energies for diamond are fairly well reproduced. However, the elastic constants as a function of temperature [Gao, et al. J. Chem. Phys., 125 (2006) 144506] and the force required to break covalent bonds were subsequently shown to be incorrectly reproduced [Pastewka, et al. Phys. Rev. B, 78 (2008) 161402(R)].

We recently developed a new bond-order potential for systems containing Si, C, and H, such as organosilicon molecules, solid Si, solid C, and alloys. This potential is based on the REBO2 potential [Brenner et al., J. Phys. C, 14 (2002) 783] for hydrocarbons and REBO for Si [Schall et al., Phys. Rev. B, 77 (2008) 115209]. Modifications to the hydrocarbon REBO potential were made to improve the description of three-atom type systems. The widespread use of Brenner’s REBO potential, its ability to model a wide range of hydrocarbon materials, and the existence of parameters for several atom types were the motivating factors for obtaining the Si-C-H (2B-SiCH) parameterization. We will briefly outline the 2B-SiCH potential and present recent simulations that examined adhesion and friction between Si and SiC tips and diamond surfaces.

The adaptive intermolecular REBO, or AIREBO, potential includes intermolecular interactions between non-bonded atoms and torsional interactions associated with a connected sequence of 3 bonds [Stuart, et al. J. Chem. Phys. 112, 6472 (2000)]. This makes the AIREBO potential capable of studying interfacial/tribological systems where chemical reactions take place. A method for extending charge transfer to BOPs, known as the bond-order potential/split-charge equilibration (BOP/SQE) method [Mikulski, et al. J. Chem. Phys. 131, 241105 (2009)], was recently integrated into a new BOP for interactions between O, C, and H. This qAIREBO potential [Knippenberg et al., J. Chem. Phys. 136, 164701 (2012)] is able to model chemical reactions where partial charges change in gas- and condensed-phase systems containing O, C, and H. The BOP/SQE method prevents the unrestricted growth of charges, without adding significant computational time, because it makes use of a quantity which is calculated as part of the underlying covalent portion of the potential, namely, the bond order. We will outline some of our recent work that uses the qAIREBO potential to examine the humidity dependence of friction in carbon-based solid lubricants.

Thermal transport in two-dimensional silicon calls for accurate force field
Ming Hu, Xiaoliang Zhang, and Kefei Li

Download talk as PDF.

Recently there has been enormous progress in the production of atomic layers that are strictly two-dimensional (2-D) and can be viewed as individual planes of atomic-scale thickness pulled out of bulk materials. Silicon is a leading material in semiconductor industries because of its wide usage in transistors, solar cells, and other electronic devices. Silicene, the silicon equivalent of graphene, has recently attracted significant attention.

In this presentation, first I will describe our recent molecular dynamics (MD) simulation on thermal transport of atomically thin silicene with improved Stillinger-Weber potential. We found that the silicene has intrinsically extremely low thermal conductivity (~ 10 W/mK at 300 K), which is about 15 times smaller than bulk Si and two orders of magnitude less than that of its carbon counterpart graphene. By quantifying the relative contribution from different phonon polarizations, we show that the phonon transport in silicene is dominated by the longitudinal modes, not the out-of-plane flexural modes as opposed to graphene.

Second, I will demonstrate the effect of surface functionalization on the thermal transport of silicene. We discovered that, contrary to its counterpart of graphene and despite the similarity of their honeycomb lattice structure, silicene exhibits an anomalous thermal behavior upon surface functionalization: the thermal conductivity of silicene significantly increases with applied hydrogen percentage.

Third, by conducting a series of non-equilibrium molecular dynamics simulations, we show that substrates have great effect on the thermal conductivity of silicene. More importantly, substrates can tune the thermal conductivity of silicene in a broad range, which is very important for the thermal management of electronic devices involving silicene. Despite the above findings, the thermal transport in two-dimensional silicon still calls for more accurate interatomic potential to precisely predict the phonon behavior in such structure. With accurate potential MD simulation can provide a guide of how to modulate the thermal transport properties of two-dimensional Si with nanoengineering and may be of use in tuning their electronic and optical properties for electronic, thermoelectric, photovoltaic, and opto-electronic applications.

Towards a physically sound description of zinc oxide (nano-)systems: SCC-DFTB – bridging a gap?
Stefan E. Huber, Matti Hellström, Michael Probst, Kersti Hermansson, and Peter Broqvist

Zinc oxide (ZnO) has a wide range of applications due to its unique catalytic, photo-catalytic, electronic and many other properties. For most of these properties, surface chemistry plays a crucial role, and thus, due to the enhanced surface-to-volume ratio, ZnO nano-systems are regarded as promising candidates for a variety of (future) technological applications such as for water-splitting, the synthesis of hydrogen peroxide or the reduction of graphene oxides to name just a few. For computational modeling, ZnO nanosystems represent a difficult task. Whereas density functional theory (DFT) methods are associated usually with severe limits in system size due to their computational demand, computationally cheaper force-field approaches lack in the description of electronic properties. In addition, the parameterization of force fields can be a tedious task. This is also related to the computational demand associated with the electronic structure methods usually used for the calculation of reference data against which force fields are then parameterized. Here, an attempt to bridge the gap between the quantum mechanical and atomistic scale, shall be presented. The goal is to use an approximate DFT method, namely density-functional tight-binding with self-consistent charges (SCC-DFTB) [1], as an intermediate step in the parameterization procedure. In particular, DFT is used for the parameterization of SCC-DFTB and the latter is then used for the major part of the parameterization of force-fields. This idea is based on the fact, that usually relatively little reference data is required to parameterize SCC-DFTB in contrast to the considerably larger training sets usually required for the parameterization of force fields. Thus, use of SCC-DFTB as an intermediate step can yield a considerable speed-up of the overall parameterization procedure. However, as an approximate DFT method, the reliability of SCC-DFTB relies per se heavily on the quality of the used parameters.

In this presentation, the focus is on the recent development, assessment and application of such an SCC-DFTB set of parameters for ZnO systems [2, 3]. It is found that SCC-DFTB in conjunction with the new set of parameters yields results in agreement with experimental and theoretical data concerning energetic, geometrical and electronic properties of ZnO bulk systems, low-index polar and non-polar surfaces and also surface-water interfaces. The method is thus capable of a reasonable description of ZnO nanosystems and can be set to use for the said purpose. In addition, its lower computational cost compared to DFT allows the treatment of considerably larger systems than is possible with DFT, while the method still retains the electrons in the model which allows the investigation of electronic properties. The latter feature makes the method also an interesting candidate in multiscale modeling frameworks in which the specific advantages of force fields and (approximate) electronic structure methods are combined to extract data at different spatial and temporal scales.

[1] Aradi, B.; Hourahine, B.; Frauenheim, Th; J. Phys. Chem. A 111 (2007) 5678.
[2] Hellström, M.; Jorner, K.; Bryngelsson, M.; Huber, S.E.; Kullgren, J.; Frauenheim, Th.; Broqvist, P.; J. Phys. Chem. C 117 (2013) 17004.
[3] Huber, S.E.; Hellström, M.; Probst, M.; Hermansson, K.; Broqvist, P.; Surf. Sci. 628 (2014) 50.

A structural fitting approach to develop effective potentials from ab initio simulations of amorphous systems
A. Carré, S. Halbert, J. Horbach, S. Ispas, C. Raynaud, O. Eisenstein, and W. Kob

The art of deriving classical potential consists in finding the correct balance between the complexity of the functional form whose implementation should be as simple as possible, the computational efficiency, and the reliabilty. In this talk, we will present a new methodology which uses ab initio data for parametrizing classical potentials suited for amorphous systems.

We have considered a simple and widespread empirical Ansatz, i.e. Born-Mayer-Coulomb, with the goal to investigate large time/size scales. As a first condition, a classical potential should be able to describe the atomic structure of the system under investigation. In order to optimize the parameters of the classical potential we have then used a cost function that involves the partial pair correlation functions extracted from ab initio simulations. The system considered is liquid silica, a liquid which represents a prototype of strong glass formers exibiting a network structure. From a methodological point of view, we will show, through the different comparisons between ab initio and classical MD, that this structural fitting procedure is a very promising approach to obtain new effective potentials suited for describing not only amorphous material but also crystalline properties.

We will present simulations carried out using this new potential and showing that CHIK reproduces with more accuracy a number of properties of molten and glass silica, than the well known BKS potential. Some tests performed on a-quartz also showed that the new potential can also be used to simulate crystalline systems with a good accuracy regarding properties such as asymmetric coordinates and elastic constants. Finally, we will discuss molecular dynamics simulations carried out to represent the structures of amorphous dehydroxylated silica surface in liquid and glass states, and show that the new potential gives a also fair description of the surface.

Systematic analysis and extension of embedded atom models
J. Jalkanen and M. H. Müser

Download talk as PDF.

Embedded-atom methods (EAM) are among the most popular classical force fields for pure metals and alloys. Despite the relative simplicity of the model, it is capable of reproducing the correct trends for a variety of physical quantities and has been found applicable to a wide range of structural problems. Of particular interest in this study is the transferability of the method between different coordination environments.

The EAM potentials are composed of a pair-wise repulsive term and an embedding term which functionally depends on an electronic density, usually formed as a linear combination of atom-centered radial charge distributions. Both the repulsion and the embedding terms can be chosen to reproduce a given equation of state exactly.

In this work, we replace each of these EAM components with different functional forms and analyze how well they describe a set of ab initio copper structures ranging from clusters to atomic wires, layer structures and bulk lattices. Since the degrees of freedom related to charge transfer are outside the scope of the conventional EAM description, our structures are chosen to consist only of symmetry equivalent atoms.

We find that the simple EAM potentials such as Gupta and Finnis-Sinclair give a satisfactory match against the ab initio data. The more complex component functionals improved the quality of the fit only by trace amounts. We find that letting the embedding function to not only depend on the density but also on the density gradients improves the fit significantly.

In a systematic expansion of the embedding energy in terms of the charge density gradients, we find that only a few terms are necessary to yield a model which reproduces the data within its accuracy limits. This model, which we call as systematically modified embedded atom method (SMEAM), overcomes the systematic errors in the EAM dissociation energies at low coordination numbers but still continues to share the difficulties with the compliance properties of the hypothetical copper diamond.

Molecular dynamics with on-the-fly machine learning of quantum mechanical forces
Zhenwei Li, James Kermode, and Alessandro De Vita

We present a molecular dynamics scheme where forces on atoms are either predicted by a machine learning (ML) scheme or, if necessary, computed with on-the-fly quantum mechanical (QM) calculations and added to the growing ML database [1]. The resulting force field is accurate and transferable, since QM accuracy can be enforced to any desired tolerance and database completeness is never required. The scheme is also efficient since the frequency of QM calls systematically decreases, eventually falling to zero in simple situations where no novel chemical processes or bonding geometries are encountered during the simulation.

We expect the method to be particularly useful for the simulation of processes where complex but recurring chemical steps are encountered (e.g. crack propagation [2]), which can be learned, while time-localised occurrences of new chemical bonding geometries (e.g. crack-impurity interactions [3] or reactions with environmental molecules [4]) cannot be ruled out, so that a fixed classical potential is not an option.

[1] Z. Li, J. R. Kermode and A. De Vita, Submitted (2014).
[2] J. R. Kermode, T. Albaret, D. Sherman, N. Bernstein, P. Gumbsch, M. C. Payne, G. Csányi, and A. De Vita, Nature 455, 1224 (2008).
[3] J. R. Kermode, L. Ben-Bashat, F. Atrash, J. J. Cilliers, D. Sherman, and A. De Vita, Nat. Commun. 4, 2441 (2013).
[4] A. Gleizer, G. Peralta, J. R. Kermode, A. De Vita, and D. Sherman, Phys. Rev. Lett. 112, 115501 (2014).

Halogen bonds, sigma-holes and molecular mechanics of modern drug candidates
Michal H. Kolar, Paolo Carloni, and Pavel Hobza

Download talk as PDF.

Pair-additive atom-centered force fields are widely used in modeling biochemical phenomena in aqueous phase such as biomolecular conformational changes or inhibitor-enzyme binding. For a long time, halogens has been playing a prominent role in drug discovery and nowadays it is estimated that as much as 35 % of approved drugs contain a halogen atom.

Halogen bond is a recently recognized type of noncovalent interaction between a halogen and a Lewis base (e.g. carbonyl oxygen). This interaction facilitates binding of numerous organic molecules (drug candidates) to proteins related to certain types of cancer, HIV infection or diabetes. Using classical force fields halogen bonding is impossible to model, because of an inherent charge anisotropy of the halogen, called σ-hole.

We aim at extending the capabilities of the current biomolecular classical force field to faithfully describe halogen bonding. An extended analysis of halogen σ-holes has been done in order to get insight into what the σ-hole actually is from spatial point of view. Using off-centered massless charge we developed a σ-hole model (called explicit σ-hole) which notably improves both the biomolecular geometries as well as energetics. The application of ESH model goes, however, beyond the biomolecular simulations, since halogen bonding is a widely used concept also in self-assembling materials and crystal engineering.

  • M. H. Kolar, P. Carloni, P. Hobza, Phys. Chem. Chem. Phys., 2014, doi: 10.1039/C4CP02621G.
  • M. Kolar, P. Hobza, A. K. Bronowska, Chem. Commun. 2013, doi: 10.1039/C2CC37584B.
  • M. Kolar, P. Hobza, J. Chem. Theory Comp. 2012, doi: 10.1021/ct2008389.

Screened empirical bond-order potential for Si-C-H
Lars Pastewka

First nearest-neighbor models are routinely used for atomistic modeling of covalent materials. Neighbors are usually determined by looking for atoms within a fixed interaction range. While these models provide a faithful description of material properties near equilibrium, the limited interaction range introduces problems in heterogeneous environments and when bond-breaking processes are of concern.

I demonstrate that the reliability of reactive bond-order potentials is improved by using an environment-dependent first nearest-neighbor definition. I will discuss general considerations and show specific results obtained using an augmented version of the second generation reactive empirical bond-order potential for C-H and the recent Si-C-H incarnation by Schall and Harrison, as well as results obtained for simple Si-C Tersoff-type potentials.

  • Pastewka, Pou, Perez, Gumbsch, Moseler, Phys. Rev. B 78, 161402(R) (2008).
  • Pastewka, Mrovec, Moseler, Gumbsch, MRS Bull. 37, 493 (2012).
  • Pastewka, Klemenz, Gumbsch, Moseler, Phys. Rev. B 87, 205410 (2013).

Quantum Chemical Topology: Towards a novel protein force field with polarisable multipolar electrostatics
P. LA. Popelier, T. Fletcher, S. J. Davie, and M. Mills

Quantum Chemical Topology (QCT) [1,2] is the collective modern name for all approaches that use the central idea of partitioning a quantum function by means of its gradient vector field. QCT creates the topological atoms of the “Quantum Theory of Atoms in Molecules” [3-5], which serves as the basis of this novel force field. Topological atoms are finite-volume, malleable boxes that do not overlap nor leave gaps between them; they exhaust space and form a mosaic of complementary shapes. The figure shows peptide-capped tryptophan with topological atoms (C: gold, N: blue, O: red, H: white) superimposed over the molecular graph, showing the intramolecular hydrogen bonds via dotted lines. This picture was generated using in-house methodology developed earlier [6,7].

We first discuss how QCT defines bonds [8]: how to draw a molecule from a molecular wave function. The quantity that is at the basis of this important realization is Vxc, a measure of interatomic exchange-correlation energy, also interpretable as a measure of covalency. This quantity, alongside with atomic self-energies (and their deformation) [9] is also crucial in a modern interpretation of stereo-electronic effects (thereby replacing concepts such as steric hindrance and (hype)conjugation) Conveniently, such energies are transferable (for 1, n interactions in saturated linear hydrocarbons) and can provide an accurate estimation of the covalent-like contribution between pairs of given interacting topological atoms A and B.

Then we focus on the electrostatic interaction, with most detail, and then on the non-electrostatic part. If atomic coordinates change then the shapes of the atoms change too, as well as their multipole moments. This complex relationship is captured by a machine learning technique called kriging. Here we show how these ideas [10] can be used to enhance the realism of the electrostatic energy [11-13], and put polarisation and charge transfer on the same footing, without having a polarisation catastrophe.

[1] Popelier, P. L. A.; Brémond, É. A. G. Int.J.Quant.Chem. 2009, 109, 2542.
[2] Popelier, P. L. A. In Structure and Bonding. Intermolecular Forces and Clusters, Ed, D.J.Wales; Springer: Heidelberg, Germany, 2005; Vol. 115, p 1.
[3] Bader, R. F. W. Atoms in Molecules. A Quantum Theory.; Oxford Univ. Press: Oxford, Great Britain, 1990.
[4] Popelier, P. L. A. Atoms in Molecules. An Introduction.; Pearson Education: London, Great Britain, 2000.
[5] Popelier, P. L. A. In The Nature of the Chemical Bond Revisited; Frenking, G., Shaik, S., Eds.; Wiley-VCH, Ch.8: 2014, p 271.
[6] Rafat, M.; Devereux, M.; Popelier, P. L. A. J. Mol. Graphics Modell. 2005, 24, 11.
[7] Rafat, M.; Popelier, P. L. A. J.Comput.Chem. 2007, 28, 2602.
[8] Garcia-Revilla, M.; Francisco, E.; Popelier , P. L. A.; Martin-Pendas, A. M. ChemPhysChem 2013, 14, 1211.
[9] Pendas, A. M.; Blanco, M. A.; Francisco, E. J.Comput.Chem. 2009, 30, 98.
[10] Popelier, P. L. A. AIP Conf.Proc. 2012, 1456, 261.
[11] Mills, M. J. L.; Popelier, P. L. A. Theor.Chem.Acc. 2012, 131, 1137.
[12] Kandathil, S. M.; Fletcher, T. L.; Yuan, Y.; Knowles, J.; Popelier, P. L. A. J.Comput.Chem. 2013, 34, 1850.
[13] Fletcher, T.; Davie, S. J.; Popelier, P. L. A. J. Chem. Theory Comput. 2014,

Building polarizable force fields for ionic materials and liquids
Mathieu Salanne

Many modeling problems in materials science involve finite temperature simulations with a realistic representation of the interatomic interactions. These problems often necessitate the use of large simulation cells or long run times, which puts them outside the range of direct first-principles simulation. In ionic systems, it is possible to introduce physically motivated model potentials for the interactions, in which additional degrees of freedom provide a ’cartoon’ of the response of the electronic structure of the ions to their changing coordination environments and allow the compact representation of many-body contributions to the interaction energy. These potentials may then be parameterized by fitting the predicted forces and multipoles to a large body of information generated from first-principles calculations. The resulting potentials are predictive, of first-principles accuracy, and have a high degree of transferability between different systems.

In this work, interaction potentials were developed for various ionic systems ranging from molten salts (fluorides and chlorides), oxides (mainly silicates and borates) and a room temperature ionic liquid (ethyl-methylimidazolium chloride mixtures with aluminium chlorides). These potentials were validated on a wide range of properties (spectroscopy experiments, activity coefficients, electrical conductivity and viscosity data), and they now allow the prediction of many other properties which remained unknown.

Polarizable force-fields for oxides (including water)
Sandro Scandolo, Paul Tangney, Nicola Seriani, Calos Pinilla, Otto Gonzalez, and Ali Hassanali

Accurate force-fields parameterized on ab initio calculations are a promising way to extend the size and time constraints of ab initio simulations. Extensive work on polarizable models for oxides has shown that accurate potentials can be constructed that reproduce faithfully the ab initio results on homogeneous bulk phases. Compounds studied so far include SiO2, MgO, MgSiO3, TiO2, and water. The force-field developed for water is all-atom, dissociable, and capable of describing with reasonable accuracy also the dissociated state.

First principles parametrized force fields for porous coordination polymers and related systems: Development and application
Rochus Schmid

Details about the development and in particular the parametrization of the force field MOF-FF [1] will be presented. MOF-FF has been designed to describe porous coordination polymers like MOFs and related systems by a bonded but flexible molecular mechanics potential. The parameters are generated from DFT reference data of model systems with the idea of "sacrificing transferability for accuracy" in mind. A Genetic Algorithm and a redundant internal coordinate based objective function are used here. In the presentation the use of such a force field for the structure prediction of framework materials as an important application is briefly sketched.

[1] S. Bureekaew, S. Amirjalayer, M. Tafipolsky, C. Spickermann, T. K. Roy, R. Schmid, Phys. Stat. Sol. B 2013, 250, 1128-1141 (DOI:10.1002/pssb.201248460).

Critical assessment of SQE model for the development of polarizable force fields
Konstantin S. Smirnov

The contribution presents results of using the split-charge equilibraion (SQE) model as a background of polarizable force field. The discussed results mainly concern two systems, silica and water, that are widely used as prototypes of semi-ionic solids and of H-bonded molecular liquids, respectively. Strengths and weaknesses of the model in the description of specific characteristics are revealed by systematic comparison with reference ab initio data and with other effective potential models.

Reactive MD simulations: From nanoscale electrochemistry to detonation initiation
Alejandro Strachan, Nicolas Onofrio, Mathew Cherukara, and Mitchel Wood

This presentation will discuss recent developments in reactive MD simulations to describe electrochemical processes and exemplify its use with two applications where chemical reactions are driven by different external stimuli.

Nanoscale resistance-switching cells that operate via the electrochemical formation and disruption of metallic filaments that bridge two electrodes are among the most promising devices for post-CMOS electronics. Despite their importance, the mechanisms that govern their remarkable properties are not fully understood, limiting our ability to assess the ultimate performance and scalability of this technology. We present the first atomistic simulations of the operation of conductive bridging cells using reactive MD with a charge equilibration method extended to describe electrochemical reactions. The simulations predict the ultrafast switching observed in these devices, with timescales ranging from hundreds of picoseconds to a few nanoseconds for devices consisting of Cu active electrodes and amorphous silica dielectrics and dimensions corresponding to their scaling limit. We find that single-atom-chain bridges often form during device operation but they are metastable with lifetimes below a nanosecond. The formation of stable filaments involves the aggregation of ions into small metallic clusters followed by a progressive chemical reduction as they become connected to the cathode. Contrary to observations in larger cells, the nanoscale conductive bridges often lack crystalline order. An atomic-level mechanistic understanding of the switching process provides new guidelines for materials optimization for such applications and the quantitative predictions over an ensemble of devices provide insight into their ultimate scaling and performance.

Detonation initiation of RDX under dynamical loading. Defects play a central role in the initiation of chemical reactions in high energy density materials under dynamical loading by spatially localizing the shock energy. Despite many decades of research the details of how the energy in the shockwave leads to self-sustained chemical reactions remain obscure. We used large-scale reactive MD simulations to identify the critical pore size that results in sustained chemistry for the nitramine RDX (C3H6N6O6) under shock loading with a piston velocity of 2 km/s. The mechanical shock causes pore collapse and the initial chemical reactions occur within picoseconds of this event and under local non-equilibrium conditions. Following these initial reactions the chemical front grows most rapidly into the collapsed pore (i.e. in the upstream direction) which we attribute to the heat lensing effect of the pore and the increased sensitivity of the amorphous material. Above the critical pore size we see the formation of final, exothermic, products that leads to significant local heating which accelerates the chemical reactions. These simulations provide atomic insight into the coupling of dynamical mechanical loads and chemical reactions that occur under extreme conditions of temperature and pressure and away from local equilibrium.

Accuracy and transferability of GAP models for tungsten
Wojciech J. Szlachta, Albert P. Bartok, and Gabor Csanyi

Download talk as PDF.

We introduce interatomic potentials for tungsten in the bcc crystal phase and its defects within the Gaussian Approximation Potential (GAP) framework, fitted to a database of first principles density functional calculations. We investigate the performance of a number of models based on a series of databases of increasing coverage in configuration space and showcase our strategy of choosing representative small unit cells to train models that predict properties only observable using thousands of atoms. The most comprehensive model is then used to calculate properties of the screw dislocation, including its structure, the Peierls barrier and the energetics of the vacancy-dislocation interaction.

All software and data are available at

OpenKIM and the future of force fields: Citability, portability and transferability
Ellad B. Tadmor

Atomistic simulations using empirical force fields play a key role in realistic scientific and industrial applications. This talk describes a US NSF-funded effort to develop an open-source online tool for promoting the use and reliability of interatomic models (force fields). The Open Knowledgebase of Interatomic Models (OpenKIM) ( allows users to compare model predictions with reference data, to generate new predictions by uploading simulation test codes, and to download models conforming to an application programming interface (API) standard that has been developed in collaboration with the atomistic simulation community. Interatomic models are archived in OpenKIM with permanent citable IDs making it possible for others to reproduce published work. Models downloaded from OpenKIM can be seamlessly used with any of a growing number of major simulation codes that support the KIM API. An overview will be given of the OpenKIM project and its main components which include the KIM API, the KIM data structure for representing arbitrary material properties, the KIM processing pipeline, the KIM visualization framework, and "RATE", a method for gauging model transferability.

Applications of ReaxFF to complex materials and material interfaces
Adri van Duin

The ReaxFF method provides a highly transferable simulation method for atomistic scale simulations on chemical reactions at the nanosecond and nanometer scale. It combines concepts of bond-order based potentials with a polarizable charge distribution. Since it initial development for hydrocarbons in 2001, we have found that this concept is transferable to applications to elements all across the periodic table, including all first row elements, metals, ceramics and ionic materials. For all these elements and associated materials we have demonstrated that ReaxFF can accurately reproduce quantum mechanics-based structures, reaction energies and reaction barriers, enabling the method to predict reaction kinetics in complicated, multi-material environments at a relatively modest computational expense.

This presentation will describe the current concepts of the ReaxFF method, the current status of the various ReaxFF codes, including parallel implementations. Also, we will present and overview of recent applications to a range of materials of increasing complexity, with applications to combustion, catalysis, aqueous phase chemistry and material failure.

How to ensure the accuracy of polarizable force fields?
Toon Verstraelen, Steven Vandenbrande, and Paul W. Ayers

In many force-field simulations, host-guest interactions determine the relevant results, e.g. the spatial distribution of guest molecules in a porous material, or the binding free energy of a ligand in an enzyme. Because the accuracy of non-covalent interactions is so critical, there is a continued interest in explicit treatments of electronic polarization in force fields since the early nineties. These so-called polarizable force fields (PFFs) seek to improve the accuracy of non-covalent interactions with a modest extra computational cost. The mutual electronic polarization of two compounds leads to an attractive force that is unattainable with few-body terms. Despite intensive research, the development of a quantitatively accurate PFF is still a very complex and labor-intensive procedure: parameters are calibrated empirically with diverse sets of reference data, including hydration energies, NMR data and ab initio data.

To facilitate the development of accurate PFFs, we propose to eliminate the empirical factors and to construct these models rigorously from an ab initio description of the electronic linear response. To this end, we recently derived “Atom-Condensed Kohn-Sham DFT approximated to Second order” (ACKS2) and showed that classical PFFs are lacking important non-local contributions due to the electronic kinetic energy. In this work, we generalize the theoretical foundations of the ACKS2 model, building on other theories than KS-DFT and including atomic multipoles to arbitrary order. We validate our approach numerically, showing that the generalized ACKS2 model reproduces ab initio reference data without any empirical calibration of parameters.


The posters will be mounted on movable walls provided by the organizers. The maximum size of a single poster should not exceed 90cm widths and 145cm height (portrait format).

The number in brackets in front of the title is the number of the movable wall where to place the poster for poster session.

[1] Ab initio study of defected surfaces in CoPt L10 alloy
Samy Brahimi and Hamid Bouzar

Download poster as PDF.

The equiatomic CoPt alloy is a promising candidate for ultra high magnetic data recording, due to its tetragonal L10 structure which has a high uniaxial magnetic anisotropy. In this work we present ab initio study of magnetic properties of defected surfaces in CoPt L10 alloy, especially surfaces with localized faults, such as antisites. Our results show that the Co (Pt) local magnetic moment is increased in the presence of Pt (Co) antisite.

Keywords: Ab initio, equiatomic alloy, CoPt L10, antisites, uniaxial magnetic anisotropy.

[2] Modeling charging and discharge of electrochemical cells from atomistic principles
Wolf Dapp

Calculating the Coulomb forces in Molecular Dynamics (MD), the effective atomic charge on each atom must be known accurately. While there are several approaches (QE, AACT) for this task, they each have shortcomings. Split Charge Equilibration (SQE) is a method combining the beneficial aspects of both QE and AACT, while avoiding their deficiencies, such as predicting a wrong charge distribution among dissociation of dimers.

To date, partial atomic charges are either assumed constant or assigned according to minimization principles that are smooth and unique functions of the instantaneous atomic coordinates. Thus, charge-transfer history effects cannot be accounted for. Such occur for instance during contact-induced charge transfer between two solid bodies. In particular, redox reactions involve a quasi-discontinuous change of the electronic state, without significant atomic rearrangement.

RedoxSQE, an extension to SQE, assigns discrete oxidation states to each atom. This enables us to model redox reactions. This talk will show proof-of-concept applications, including a complete atomistic nano battery, which self-consistently shows many generic properties of macroscopic batteries. Among other applications, redoxSQE can help to better understand the processes happening at the electrode-electrolyte interface of a real battery.

[3] Molecular dynamics simulations of edge dislocations in iron interacting with grain boundaries
A. Gubbels-Elzas and B. J. Thijsse

Download poster as PDF.

Plasticity in iron is governed by the motion of dislocations. Grain boundaries in the material, whether between similar or between dissimilar grains, act as a barrier for dislocation motion. Dislocations can be reflected by a grain boundary, absorbed, the direction of motion can be changed or a combination of these phenomena can occur. The influence of a grain boundary between a bcc and a fcc grain on the motion of an edge dislocation is studied with molecular dynamics simulations. The bcc material is described by an iron-potential capable of describing dislocation behaviour (L. Malerba et al., Journal of Nuclear Materials, 406, 2010). To study the influence of a grain boundary between a bcc and a fcc grain, this iron-potential is combined with a nickel-potential (G. Bonny et al., MSMSE, 17, 2009) to describe both phases. The interaction of edge dislocations with differently oriented grain boundaries is studied.

[4] Modeling the effect of occupational disorder in Lithium-Titanium-Oxide materials
Hendrik Heenen, Saskia Stegmaier, Christoph Scheurer, and Karsten Reuter

Li4Ti5O12 lithium-titanium-oxide (LTO) materials represent an interesting alternative to conventional graphite based anodes for Li ion batteries, due to their high cycling stability, safe operation at high working potentials and a fast charge-discharge behavior [1,2]. In order to rationalize and improve the material's properties and performance, a detailed understanding of the underlying atomic scale features and processes is required. In this context, different studies have been conducted to evaluate transition state energies and diffusion pathways in LTO [3-6]. The structural models of LTO used in these studies are all lacking a proper treatment of the statistics of the system's configuration space. The spinel structure phase of Li4Ti5O12 exhibits mixed occupancy of octahedral sites by Li and Ti ions in the LTO bulk phase which allows for a high degree of occupational disorder. In order to sample a thermodynamically valid representation of said disorder, it is necessary to consider larger structural supercells.

In the present study force field calculations have been used and validated by DFT data to get a glimpse of the system's vast configuration space and elucidate the impact of its proper sampling. This approach will be extended to a full multiscale representation of LTO in the future.

[1] M.V. Reddy, G.V. Subba Rao, B.V.R. Chowdari, Chem. Rev., 2013, 113(7), 5364-457.
[2] M.K. Song, S. Park, F.M. Alamgir, J. Cho, M. Liu, Mat. Sci. Eng. R, 2011, 72(11), 203-252.
[3] Z. Zhong, C. Ouyang, S. Shi, M. Lei, Chem. Phys. Chem., 2008, 9(14), 2104-8.
[4] H. Shiiba, M. Nakayama, M. Nogami, Solid State Ionics, 2010, 181(21-22), 994-1001.
[5] M. Vijayakumar, S. Kerisit, K.M. Rosso, S.D. Burton, J.A. Sears, Z. Yang, G.L. Graff, J. Liu, J. Hu, J. Power Sources, 2011, 196(4), 2211-2220.
[6] B. Ziebarth, M. Klinsmann, T. Eckl, C. Elsässer, Phys. Rev. B, 2014, 89, 174301.

[6] Systematic modifications of embedded-atom potentials
J. Jalkanen and M. H. Müser

A good transferability of a force field between environments of varying coordination numbers is generally a difficult goal to achieve. It is commonly said that an excellent description of bulk quantities of a given material requires sacrifices in terms of surface properties, and vice versa. The problem is complicated by the fact that often a change in the atomic neighborhood tends to be accompanied by a redistribution of charge, leading to a degradation in the quality of a given classical model potential for the material.

To study the applicability of embedded-atom method-type potentials (EAM) in this situation, we construct a set of copper structures where all charge transfer is disallowed by symmetry. Using ab initio calculations for this set as our reference data, we vary the shapes of the terms constituting a generic EAM potential to identify the functions with the best transferability properties.

We find that most choices for the constituent functions lead to a similar quality of description. However, a significant improvement is observed when the embedding energy is allowed to depend on the gradients of the charge density. Using the Gupta potential as our reference, we perform a systematic test of including the lowest order combinations of charge density gradients to the embedding energy functional. According to our results, only a few terms are needed to improve the fit to a level where it is on a par with the accuracy of the ab initio results, without a significant computational overhead.

We show that this force field, which we refer to as systematically modified embedded atom method (SMEAM), improves trends of the nearest neighbor distances and dissociation energies of the structures as a function of the coordination. Additionally, SMEAM reproduces the equations of state for several bulk structures with an excellent accuracy, and gives also reasonable values for the fcc surface- and vacancy energies.

[8] SQE-based polarizable water model implemented in GROMACS
Yuriy Khalak and Mikko Karttunen

Split Charge Equlibration was implemented into GROMACS, a popular molecular mechanics toolkit, to allow for easy modeling of polarization in biological soft-matter systems. Water was chosen as a proof of concept system and was parametrized based on the TIP4P/2005 model. Shielding effects were included by approximating the overlap integral of Slater orbitals with a sigmoid function to reduce the computational load.

[9] QM/MM assessment of the proton and electron affinities of ligated heme iron centers in various enzyme families
Anikó Lábas, Balázs Krámos, Tibor Szilvási, and Julianna Oláh

Heme, an iron porphyrin, is one of the most important prosthetic groups in metalloproteins. It is the active agent of various enzyme families, e.g. if cytochrome P450s, peroxidases and globins. These proteins show striking diversity in their mode of action, and are involved in many biological processes, like drug metabolism or hormone synthesis. A remarkable feature of cytochrome P450 enzymes is their ability to oxidize inert carbon-hydrogen bonds, which keeps on inspiring the design of novel catalysts. The generally accepted mechanism of cytochrome P450s is called “rebound mechanism”: compound I abstracts a hydrogen from the substrate to form the hydroxo complex and a substrate radical, which are rapidly recombined to generate the hydroxylated product.[1] The ability of the FeIV-oxo complex to abstract a hydrogen is facilitated by the electron push effect of the proximal ligand and is related to the bond energy of the forming O-H bond (D(O-H)) of the hydroxo complex. This bond energy can be approximated by an empirical form,[2] containing the pKa of hydroxo complex and the one-electron reduction potential of compound I: D(O-H) = 23.06 • E0compound I + 1.37 • pKa,hydroxo complex + constant

The goal of this computational study is to compare the relationship between the proton and electron affinities of the heme centres of various enzyme families (CYPs, APOs, HRP, LiP) and to examine the electron donating effect of the different proximal ligands. For this, molecular dynamics simulations were carried out in the hydroxo complex state, and QM/MM geometry optimizations of a large number of protein conformations were carried out in order to collect statistically significant data that can be used to estimate the pKa and redox potential values of the various enzymes.

Acknowledgement: This work was supported by Gedeon Richter Plc. and OTKA Grant No 108721.

[1] John T. Groves, Nat. Chem.2014, 6, 89–91.
[2] James M. Mayer, Acc. Chem. Res. 1998, 31, 441-450.

[10] Open collaboration that uses NMR data to judge the correctness of phospholipid glycerol and head group structures in Molecular Dynamics simulations
Patrick F. J. Fuchs, Matti Javanainen, Antti Lamberg , Markus S. Miettinen, Luca Monticelli, Jukka Määttä, O. H. Samuli Ollila, Marius Retegan, and Hubert Santuz

We compare the C-H order parameters measured by Nuclear Magnetic Resonance (NMR) experiments to those predicted by 12 different molecular dynamics (MD) simulation models. We focus on the order parameters of the lipid headgroups and glycerol backbones in phospholipid bilayers.

Only two of the models (CHARMM36 [1] and Maciejewski-Rog [2]) give a reasonable agreement with experiments for a fully hydrated lipid bilayer.

We then compare (for the two best-performing models at full hydration and for the Berger model [3], the most used lipid model in the literature) to NMR experiments the changes in the order parameters as a function of hydration level, NaCl and CaCl2 concentrations, and cholesterol content. The results clearly show that the glycerol and headgroup structures in the Berger model are not realistic, the Na ion partitioning is significantly too strong and cholesterol-induced structural changes are overestimated. The CHARMM36 and Maciejewski-Rog perform better, but the Na partitioning is too strong at least in the latter.

This is an open science project that is progressed at All the results and discussions are available at that address.

[1] J. B. Klauda, ..., R. W. Pastor. J. Phys. Chem. B 114 7830 (2010).
[2] A. Maciejewski, M. Pasenkiewicz-Gierula, O. Cramariuc, I. Vattulainen, T. Rog. J. Phys. Chem. B 118 4571 (2014).
[3] O. Berger, O. Edholm, F. Ja¨hnig. Biophys. J. 72 2002 (1997).

[12] The structure and dynamics of the interface between the edges of clay nano-particles and aqueous solutions by molecular modelling
Maxime Pouvreau and Andrey Kalinichev

Clay minerals consist of layered tetrahedral silicate sheets and octahedral hydroxide sheets, where the octahedral sites are occupied by Al in the unsubstituted structure. In many cases, various so-called isomorphic substitutions are possible on both tetrahedral (e.g., Al3+ for Si4+) for and octahedral (e.g., Mg2+ for Al3+; Fe2+ for Al3+, etc.) sites. This leads to a very significant compositional and structural diversity of clays. The substitutions create a negative charge on the clay layer, which is usually compensated by various interlayer cations (alkali metals, alkaline earths, etc.) which can be easily hydrated. Because of the potential large-scale use of clay minerals as potential host rocks for radioactive waste repository, good understanding of the molecular mechanisms controlling radionuclide sorption and mobility in clays is necessary.

Clay surfaces can be divided into basal (parallel to the layer), and edge surfaces. When the cleavage of the clay crystal is imperfect, the sites are heterogeneously protonated in the presence of water. The figure illustrates a typical case of the edge surfaces in pyrophillite, the unsubstituted clay structure. In this way, as an addition to the ability of all surfaces to complex with various aqueous species, edge surfaces are also the place of acid-base reactions.

While molecular simulations involving classical force fields can be used to study chemically inert surfaces, quantum chemistry calculations are necessary to fully describe their reactivity. In the literature, such calculations deal with the acidity of edge sites, and with the structure and dynamics of water and complexes at the edges.

Using ab initio molecular dynamics simulations and geometry optimizations, the goals of this PhD work are to improve the classical force field ClayFF parameters related to the surface hydroxyls and to further describe the molecular processes (acid-base reactivity, complexation) occurring at clay edges.

[13] On the comparison of different barostat implementations for the prediction of the breathing behavior in MIL-53 frameworks
Sven Rogge, Louis Vanduyfhuys, Toon Verstraelen, Guillaume Maurin, and Veronique Van Speybroeck

Metal-organic frameworks (MOFs) are microporous crystalline materials, compromising metallic clusters or chains connected by organic linkers, forming a scaffold-like network with cavities and/or channels, and were first synthesized in 1999 [1]. Some of these frameworks, such as MIL-53, may undergo reversible, large-amplitude structural deformations. These transitions can be induced by applying a mechanical pressure [2], paving the way for mechanical energy storage applications.

However, these applications require a clear understanding of the response of the material on the applied pressure, especially the pressure at which the transition is induced. We determine this transition pressure with the aid of Molecular Dynamics (MD) simulations in the NPT ensemble, allowing for anisotropic volume transformations. Since the time needed for these materials to undergo these large-amplitude transformations is on the order of several hundreds of picoseconds, a complete quantummechanical treatment of the forces is not yet feasible. Therefore, force fields, specifically aimed for these MOFs and developed at the Center for Molecular Modeling [3], are used in conjunction with a simulation package that focuses on using these force fields [4].

To model the NPT ensemble, we have implemented the Nosé-Hoover chain thermostat [5-7], while the barostat implementation of Martyna, Tobias and Klein is used [8]. The simulation results are compared to the Berendsen [9] and Langevin barostats [10], as well as to experimental results from our partners in Montpellier [11].

From these MD simulations, it is observed that pressure fluctuations are significant, and may be up to four orders of magnitude larger than the average applied pressure. As a consequence, one may expect that transitions can be driven by these fluctuations, rather than by the mean applied pressure. We propose a method to distinguish between both possibilities, allowing to successfully determine the transition pressure. Moreover, the effect of the Berendsen and Langevin barostat on these pressure fluctuations is quantified for the MIL-53 type frameworks.

[1] H. Li, M. Eddaoudi, et al. Nature 402, 276 (1999).
[2] I. Beurroies, M. Boulhout, et al., Angew. Chem. Int. Edit., 49, 7526 (2010).
[3] L. Vanduyfhuys, T. Verstraelen, et al., J. Chem. Theory Comput., 8, 3217 (2012).
[4] T. Verstraelen, L. Vanduyfhuys, et al.,
[5] S. Nosé, J. Chem. Phys., 81, 511 (1984).
[6] W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
[7] G. J. Martyna, M. L. Klein, et al., J. Chem. Phys., 97, 2635 (1992).
[8] G. J. Martyna, D. J. Tobias, et al., J. Chem. Phys., 101, 4177 (1994).
[9] H. J. C. Berendsen, J. P. M. Postma, et al., J. Chem. Phys., 81, 3684 (1984).
[10] G. Bussi, M. Parrinello, Phys. Rev. E, 75, 056707 (2007).
[11] P. G. Yot, Z. Boudene, et al., Chem. Commun., 50, 9462 (2014).

[14] Cauchy violation and charge transfer potentials
Sergey Sukhomlinov and Martin H. Müser

In a crystal where particles interact through central pair forces, the elastic tensor Caßγξ possesses complete symmetry in its four indices, thus providing relations between components. These equalities are known as Cauchy relations. For a crystal of cubic symmetry only one equation remains

C12 = C44

In the presence of many-body interactions Cauchy relation no longer holds, and violation C44 - C12 from it should be considered instead. As a consequence, Cauchy violation could be considered as an estimate of how well many-body effects are described by a potential.

Chemical potential equalization (CPE) models [1, 2, 3] provide a promising way of constructing force fields with environment-dependent parameters. Assuming the specific CPE model is appropriate and its parameters are reliable, the force field can then be used to study a system’s behaviour at the microscopic level. In such polarizable force fields, the atomic charges serve as a linkage between the electronic distribution in the system and its instantaneous structure at the atomic level. Present contribution addresses the question of how violation from Cauchy relation and its slope in the region of small temperatures and ambient pressures is reproduced by the CPE models in ionic compounds.

As a model system we considered rocksalt, which in many cases can be solved analytically. We analyzed how two different CPE models reproduce many-body effects in NaCl structure by considering Cauchy violation. Charge equilibration models that were used in this study include electronegativity equalization model (EEM) [2] and split-charge equilibration (SQE) [3] model. It was found that for rocksalt structure compounds the functional form of Cauchy violation and of its pressure derivative is the same for all CPE models considered given the fact that bond hardness κs, which is a parameter of SQE model, does not depend on the environment. Both Cauchy violation and its pressure derivative are predicted to have strictly positive values, which agrees with the results of density functional theory computations and available experimental data. Nevertheless, it should be noted that the value of Cauchy violation was found to depend crucially on screening effects.

[1] D. M. York, W. J. Yang J. Chem. Phys. 104, 159-172 (1996).
[2] W. Mortier, S. Ghosh, S. Shankar J. Am. Chem. Soc. 108, 4315-4320 (1986).
[3] R. A. Nistor, J. G. Polihronov, M. H. Müser, N. J. Mosey, J. Chem. Phys. 125, 094108 (2006).

[15] Complete non-bonding force field derived from monomer electron densities
S. Vandenbrande, T. Verstraelen, L. Vanduyfhuys, M. Waroquier, and V. Van Speybroeck

When atomistic force-field simulations involving new porous materials or new ligands are carried out, carefully tuned nonbonding force-field parameters are not readily available. To avoid the expensive and daunting refinement of such parameters, it is common to take them from relatively old models, such as UFF or MM3. Recently, it was shown how one can refine such parameters using SAPT computations on a set of model dimers. However, SAPT computations are still expensive and not simply applicable to large molecular systems.

In this paper, we propose an efficient protocol to derive reliable nonbonding force fields from DFT computations on isolated monomers. The first step consists of an AIM analysis with the minimal-basis iterative stockholder method, which also yields a model of the valence electron density in terms of atom-centered Slater functions. In the second step, a van der Waals potential (for the Pauli, dispersion and induction interactions) is extracted from the partitioning data. The Slater functions and atomic charges are used directly to model the electrostatic interactions, including the penetration effect. The proposed method was tested on the S66x8 database of intermolecular interaction energies. Our force field shows an RMSD of less than 1kcal/mol when compared to CCSD(T)/CBS reference energies.

[16] QuickFF: Toward a generally applicable methodology to quickly derive accurate force fields for Metal Organic Frameworks from ab initio input
Louis Vanduyfhuys, Steven Vandenbrande, Toon Verstraelen, Rochus Schmid, Michel Waroquier, and Veronique Van Speybroeck

Metal Organic Frameworks (MOFs) are a very promising class of materials for various applications such as gas separation, detection and storage, shock absorbers, chemical catalysis and more. However, much research is still required before we can tap into the full potential of these materials, especially in the field of theoretical calculations. Although a lot has already been done, there is still much room for advancement.

One of the main reasons, is that these MOF materials are computationally very challenging to simulate. Full ab initio calculations are very demanding and can only describe a limited time and length scale. Therefore, it is of crucial importance to have accurate force fields at hand to extend the theoretical simulations to larger time and length scales. Most of the force fields that are currently at hand are either not accurate enough or limited to one particular system. Therefore, in this work, we do not propose a new force field itself, but a methodology to quickly derive these force fields from ab initio calculations. By doing so, we combine the accuracy of system-tailored force fields with the general applicability of universal force fields. Although this methodology, called QuickFF, can in principle be applied to any system of choice, it was mainly developed to construct force fields for MOFs. Therefore, we will illustrate its performance here mainly for MOFs. The required input data consists of an ab initio calculation of the geometry and Hessian in equilibrium and if necessary a set of non-bonding parameters (charges and van der Waals parameters). The mathematical expression for the covalent energy is kept simple to ensure robustness and to avoid fitting deficiencies as much as possible. The proposed methodology is split into three steps: the first two steps address the strong correlations between parameters while the last step is a refinement of the model for higher accuracy. The resulting force fields are shown to accurately reproduce both the ab initio geometry and frequencies in equilibrium for a large set of small to medium-sized organic molecules and outperform well-known general force fields such as UFF and GAFF. Furthermore, it is illustrated that QuickFF can be used to easily derive accurate force fields for more complex systems such as MOFs.

The methodology is implemented in a user-friendly Python code, which requires a minimum of input from ab initio calculations. As a result, accurate force fields for isolated molecules can easily be derived from ab initio calculations with only a minimal effort.

[17] Dynamics at a Janus interface
Michael von Domaros, Dusan Bratko, Barbara Kirchner and Alenka Luzar

Electric field effects on water interfacial properties abound, ranging from electrochemical cells to nanofluidic devices to membrane ion channels. On the nanoscale, spontaneous orientational polarization of water couples with field alignment, resulting in an asymmetric wetting behavior of opposing surfaces – a field-induced analog of a chemically generated Janus interface. Using atomistic simulations, we uncover a new and significant field polarity (sign) dependence of polarization dynamics in the hydration layer. Applying electric fields across a nanoparticle, or a nanopore, can lead to close to two orders of magnitude difference in response times of water polarization at opposite surfaces. Typical time scales are within the O(0.1) to O(10) picosecond regime. Temporal response to the field change also reveals strong coupling between local polarization and interfacial density relaxations, leading to a nonexponential and in some cases, nonmonotonic response.

This work highlights the surprisingly strong asymmetry between reorientational dynamics at surfaces with incoming and outgoing fields, which is even more pronounced than the asymmetry in static properties of a field-induced Janus interface.