# List of Abstracts

# Talks

**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 http://potfit.sourceforge.net. 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ömDownload 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 LiDownload 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 BroqvistZinc 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. KobThe 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üserDownload 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 VitaWe 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 HobzaDownload 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. MillsQuantum 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, dx.doi.org/10.1021/ct500416k.

**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 HassanaliAccurate 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 WoodThis 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 CsanyiDownload 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 www.libatoms.org.

**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) (https://openkim.org) 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. AyersIn 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.

# Posters

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 BouzarDownload 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. ThijsseDownload 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 ReuterLi4Ti5O12 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üserA 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 KarttunenSplit 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áhHeme, 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 SantuzWe 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 nmrlipids.blogspot.fi. 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 KalinichevClay 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 SpeybroeckMetal-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., http://molmod.ugent.be/software/

[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üserIn a crystal where particles interact through central pair forces, the elastic tensor
C_{aßγξ} 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

C_{12} = C_{44}

In the presence of many-body interactions Cauchy relation no longer holds, and violation
C_{44} - C_{12} 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 SpeybroeckWhen 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 SpeybroeckMetal 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 LuzarElectric 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.