Multi-State VALBOND for Atomistic Simulations of Hypervalent Molecules, Metal Complexes, and Reactions.

The implementation, validation, and application of the multi-state VALBOND method for transition-metal-containing and hypervalent molecules are presented. This approach is particularly suited for molecules with unusual shapes and systems that need to be described by a superposition of resonance structures, each of which satisfies the octet rule. The implementation is based on the original VALBOND force field and allows us to smoothly switch between resonance structures, each of which can be characterized by its own force field, including varying charge distributions and coupling terms between the states. The implementation conserves total energy for simulations in the gas phase and in solution and is applied to a number of topical systems. For the small hypervalent molecule ClF3, the barrier for pseudorotation is found to be 4.3 kcal/mol, which compares favorably with the experimentally measured value of 4.8 kcal/mol. A transition-metal-containing complex, cisplatin, is characterized by six resonance states, for which the vibrational spectrum is found to be in good agreement with experiment. Finally, umbrella sampling simulations of the SN2 reaction BrMe + Cl- → Br- + MeCl in solution yield a barrier height of 24.6 kcal/mol, in good agreement with experiment (24.7 kcal/mol).


Introduction
Atomistic simulations are a powerful means to investigate the structure and dynamics of complex systems, such as proteins in solution. Over the past 10 to 15 years the quality of such simulations has continuously improved not least due to the availability of increasingly accurate but empirical energy functions. Initially, their development focused on the chemical structure and dynamics of macromolecules, including peptides and proteins. [1][2][3][4][5][6][7] More recently, considerable progress has been made by incorporating effects due to anisotropic charge distributions (multipolar electrostatics) [8][9][10][11][12][13][14][15] and polarizability. [16][17][18] Such approaches have been very successful for applications in thermodynamics, computational spectroscopy or for investigating ligand binding to proteins.
However, when the complexity of the electronic structure of the compounds considered increases, as is the case for halogenated molecules or systems involving transition metals (TMs), even such more refined energy functions have their limitations. The situation is particularly demanding for metal-containing systems. Ideally, it would be possible to investigate the dynamics of such systems using mixed quantum mechanical/molecular mechanics molecular dynamics simulations (QM/MM MD). [19][20][21] In QM/MM the total system is divided into a (typically small) region for which the energy is calculated quantum mechanically -i.e. from electronic structure methods -and an environment which is treated with a conventional force field. However, due to the N 3 scaling of the secular determinants that need to be diagonalized, where N is the number of basis functions used in the QM part, routine applications can only be envisaged for systems including several tens of heavy atoms in the quantum region (fewer if metals are involved), and by using density functional theory (DFT) to describe the energetics. 22 Typical simulation times that can be afforded under such circumstances are on the order of several 10 picoseconds which is usually too short for reactions with appreciable activation energies. 23 An exception to this are semiempirical quantum methods, such as density functional tight binding (DFTB) 24,25 Hence, if the energetics and dynamics of reactive processes involving metals needs to be characterized, alternative methods to describe the energetics of the systems are required.
Transition metals play important roles for processes including electron transport, oxidation reactions, or oxygen transport. [26][27][28][29] Also, metal complexes find broad applications for both homogeneous and heterogeneous catalysis in organic synthesis. 30 However, atomistic simulations using empirical force fields for metals in particular and hypervalent compounds in general pose particular challenges. The difficulty arises due to variable coordination number (as exemplified for Co(CO) 4 , Cr(CO) 5 , or Cr(CO) 6 ), indistinct topology (π ligands can bind to metal in different ways), the different oxidation states (eg. Cu(I) and Cu(II) shows tetrahedral and square planar geometries respectively) 25,[31][32][33] and the coupling of electronic and vibrational degrees of freedom (i.e. Jahn-Teller distortion).
In the past, several methods reminiscent of empirical force fields have been developed to treat transition metal complexes and organometallics. The main idea underlying metal force fields is to extend a conventional force field with one or several additional terms to capture specific properties of TMs. These include large angle distortions including bending potentials with several minima, the effects due to d−electrons (such as the Jahn-Teller effect) or binding of TMs to the π−electron cloud (as in ferrocenes). As an example, the difference in bonding of axial versus equatorial ligands to a TM can be described by two different parameter sets for the two types of ligands. On the other hand, such an effect can be captured within the framework of ligand field stabilization as it has been done in ligand field molecular mechanics (LFMM). 34 Such an approach is also capable of capturing stereo-electronic effects of partially filled d-orbitals in coordination compounds.
Alternatively, the VALBOND force field [35][36][37][38] field, which is based on valence bond theory and hybrid orbital strength functions, is able to describe the angular distortions for a wide range around the metal center. According to the VALBOND formalism, 35-38 hypervalent molecules are molecules that exceed an octet electron count (e.g. ClF 3 or SF 6 ). Hypervalent molecules i.e. molecules where certain atoms have more than the formal 8 electrons, are treated by using resonance structures that formally fulfill the octet rule (or duodectet rule). For metal complexes which primarily use the ns and (n-1)d metal orbitals in forming bonds, electron counts in excess of 12 electrons are considered hypervalent. 39 Therefore, most common transition metal complexes are considered to be hypervalent and approaches such as VALBOND or VALBOND-TRANS are useful to follow the structural dynamics of TM-complexes. 33,37,38,[40][41][42] In the original VALBOND formalism hypervalent molecules were represented by combinations of 2-center-2-electron and 3-center-4-electron bonds; in the localized framework 3-center-4-electrons comprise resonance among two primary Lewis-like configurations. This approach was successful for hypervalent first row molecules (such as ClF 3 ) as well as TM-complexes. [35][36][37][38] Originally, it was only an angle potential. Later it was extended to model the trans-influence in metal complexes by introducing two energy and distance dependent penalty functions to the original VALBOND and the method called VALBOND-TRANS. 40 Besides the structural characterization of metal complexes and the dynamics around them, studying chemical reactions is fundamental in chemistry and biology. Since chemical reactions involve the breaking and making of chemical bonds, conventional force fields can not be used for such applications. However, over the years several methods such as reactive force field (ReaxFF), 43 empirical valence bond (EVB) 44 and its multi-state variant MS-EVB which was primarily used to investigate proton transport, 45,46 adiabatic reactive molecular dynamics (ARMD) 47,48 and multi-state adiabatic reactive molecular dynamics (MS-ARMD) 49 have been developed to follow chemical reactions in the gas phase and in solution. Here, we set out to generalize VALBOND 35-38 by developing a comprehensive framework approach based on resonance among structures that fulfill the octet rule. This new approach treats angular degrees of freedom within a valence bond description and further includes chemical bonds, dihedrals and van der Waals contributions by flexible parametrizations of individual resonance structures ("states").
The new method, referred to as multi-state VALBOND (MS-VALBOND), as discussed in the present manuscript explicitly treats differences in charge distributions between the resonance structures/states and allows bond breaking and bond formation. It is therefore a method that allows to investigate chemical reactions involving hypervalent and transition metals at approximately the speed of empirical force fields. The method was implemented in CHARMM 5 and was tested by applying it to three different model cases including the hypervalent molecule ClF 3 , a prototypical metal complex (cisplatin) and the chemical reaction of BrMe + Cl − → Br − + MeCl.

Computational Methods
Two models of bonding in hypervalent molecules have been proposed previously: the original Pauling approach expanded the valency by using extravalent orbitals (e.g., d-orbitals for main group elements) to accommodate electrons that exceed the octet. Alternatively, hypervalency has been described using an ionic resonance model in which linear arrangements of three atoms undergo three-center/four-electron (3c/4e) bonding. The latter model, which is supported by modern analyses of electronic structure, 50 55 In the NBO framework, the 3c/4e bond corresponds to a donor-acceptor interaction comprising donation of a lone pair of electrons from Y into the X-Y antibond.
The present approach addresses several limitations of empirical force fields for hypervalent molecules in general and metal-containing complexes in particular. First, the angular terms AXB (where A and B are atoms of ligand and X is a metal or hypervalent center) are treated by VALBOND which enables description of geometries for the bending degrees of freedom around metal centers. Secondly, each relevant resonance structure is treated as a different connectivity ("state"). The total energy of the system is described by the Hamiltonian where the diagonal elements are the force fields of each connectivity/resonance structure H ii and the off-diagonal elements are the couplings H ij between them. This formulation is akin to an the empirical valence bond (EVB) ansatz. 44,56 For a single molecule this offers the opportunity to treat the same chemical bond with different parametrizations so that effects such as alteration between single, double or partial-double bonds of one and the same bond depending on the environment can be captured. This information would be available in the strength (force constant k) and equilibrium bond length r e . It should also be possible to represent different oxidation states, e.g. Cu(I) and Cu(II). 25 Furthermore, the transinfluence 40 (known from metal-complexes) which was incorporated in VALBOND-TRANS, can be readily included. This can be achieved through suitable parametrization of the offdiagonal elements which is described further below.
Diagonalization of H yields n eigenvalues and eigenvectors and the lowest eigenvalue E 0 with is used for propagation of the system.
The contribution of state i towards the ground state energy E 0 is p i = ν i 0 ν i 0 . As all elements in H depend on internal coordinates, transitions between states along a trajectory can be followed by considering the populations p i of the states.

Parametrization of the Diagonal Elements
Each diagonal term describes a resonance structure. Thus each diagonal matrix element is treated and parametrized independently (see section 3.1). The diagonal element (H i i (X)) represents the energy of the i th resonance structure for a given configuration X. Bonds which break and form or which change dissociation energy or equilibrium structure between the states are described by Morse bonds (see equation 3).
In equation 3, r is the separation between the atoms in the bond, r e is the equilibrium distance, D e is the dissociation energy of the bond, and β describes the steepness of the short range part of the potential. Angles are described by the VALBOND formalism 35 which is rooted in valence bond theory. 57 Similarly, following previous work, 58 intramolecular van der Waals interactions are also described by Morse potentials but with small dissociation energy. This functional form has been shown to provide a better representation of the energy, especially at the short distances on the repulsive side of the energy function. Resonance structures for 3-center-4-electron bonds nearly always involve short separation in which distances (e.g., the F D − Cl A distance in resonance structure I in Figure 1) are on the repulsive side of the non-bond potential.
The non-bonded parameters of individual atoms were obtained by fitting Morse potential to the ab initio potential energy curve obtained along the Cl-He and F-He bond. Helium was used as a probe because it interacts negligibly with its environment and provides a meaningful starting point for further refinement of the parameters (see Figure S1). To describe the non-bonded interactions between non-identical atoms combination rules akin to the Lorentz-Berthelot rules are used: and r e ij = r e i +r e j 2 . It would also be possible to represent the van der Waals interactions with a conventional 12-6 form or an exponential decay with a long-range dispersion term.
To illustrate the application of the formalism the example of ClF 3 is considered, see Figure   1. The number of bonds and resonance structures required to describe the system depends on the number of valence electrons. To satisfy the octet rule one must only consider covalent Figure 1: The three resonance structures, I to III, for hypervalent ClF 3 such that each structure fulfills the octet rule. In resonance structure I, two fluorine atoms (F B and F C ) form a bond to the Chlorine atom Cl A to give octet-consistent ClF + 2 . The third atom, F D , is the octet-consistent fluoride anion F − . At the equilibrium T-shaped geometry, the F-Cl-F bond angles of the dominant resonance structures are close to 90 • , corresponding to high p-character in the Cl-F bonds of the ClF + 2 fragments.
The fully electron saturated (ionic) state formed by F − and Cl − would contain 32 electrons, hence 4 electrons are involved in bonding. Two electrons per covalent bond leads to N = 2 bonds per resonance structure. As a total of K = 3 bonds can be accommodated in ClF 3 , there are N K = 2 3 = 3 possible states (resonance structures) that need to be described, see Figure 1.
The force field description for resonance structure I consists of two bonds, which are represented each as a Morse potential with a set of fitted parameters and a VALBOND angle energy.
The angle α BAC between these bonds with the atom A at its apex is treated using valence bond theory. Within valence bond theory, for central atom A with sp m d n hybridization the energy depends only on the force constant k. 57,59,60 The general expression for the strength of hybrid orbitals S as derived by Pauling 60 is where S max is the maximum of the strength function of hybrid orbitals and ∆ is the overlap integral and their functional forms are given by In VALBOND, the overlap (∆) between two hybrid (bond-)orbitals (ψ 1 , ψ 2 ) is expressed as a function of the hybridizations m and n which are the formal p and d hybridizations. The force constant k can be chosen to best describe the bending frequency.
The nonbonded interactions in ClF 3 for resonance structure I (see with cyclical permutations for the other two states II and III. The parameters for these diagonal matrix elements (2c-2e bonds) are fitted for the resonance structure, i.e. they are fitted for ClF + 2 in the case of ClF 3 .

Parametrization of the Off-Diagonal Elements
The off-diagonal terms describe the coupling between states i and j. Within the present context, these are stabilizations between two resonance structures. The resonance between structures i and j are modeled as a scaled Morse potential depending on the angle between the bonds that are modified upon resonance from one to another resonance structure (see Figure 2) according to The reverse resonance (j → i) is modeled along the same line In order to maintain the symmetry of the Hamiltonian, the final off-diagonal element is The off-diagonal Morse parameters are then fitted to reproduce reference electronic structure data, vide infra. Effects such as the trans influence 61 or Jahn-Teller distortion 62 (coupling of vibrational and electronic degrees of freedom) can be readily included by suitably parametrizing the off-diagonal elements that reflect these effects. For example, resonance structure I is transformed by donation of a lone pair on the anion F D into the (Cl A F C ) anti-bond; such donor-acceptor interactions are the equivalent of "arrow-pushing" depictions of resonance. Shaded and empty circles indicate the lobes of the anti-bonding orbital and, for clarity, the molecular geometry shown is slightly distorted from the equilibrium T-shape. The strength of the donor-acceptor interaction depends on the overlap of donor and acceptor orbitals, hence, on the molecular geometry. At the geometry depicted above, the I→II overlap is substantially larger than the II→I overlap due to the more co-linear arrangement of the lone pair hybrid and the anti bonding orbital.
For ClF 3 the resonance between states i and j is modelled as the donation of the free electron pair into the anti-bonding orbital of the bond that is broken. 39 For the resonance between states I and II, the AC bond is broken while bond AD is formed. Considering p orbitals, the donation of the free electron pair from F D into the anti-bond depends on the angle between the three atoms A, C, and D. Scaling of the off-diagonal term (see Figure 3) is achieved with the term (cos(α DAC ) − 1) 2 which switches between a minimal amplitude (ω CAD = 0 • ) and a maximum amplitude ω CAD = 180 • , thus modeling the expected donor-acceptor overlaps.

Solute-Solvent Electrostatic Interaction in MS-VALBOND
Because including all solvent molecules into the description of the states would be computationally demanding, it is preferable to perform a mixed MS-VALBOND/MM simulation in the spirit of a mixed quantum mechanical/molecular mechanics (QM/MM) approach.
Since for each resonant structure included in the simulation the charge distribution typically differs, electron flow (or charge reorganization) needs to be followed during transitions between states. The different charge distributions characterizing each individual state lead to different solvent-solute interactions for each state and the forces need to be determined such that total energy is conserved. This concerns primarily the electrostatic interactions E elstat = Q solv q k |R−r| where the Q solv are the static solvent charges, q k are the solute charges of atom k which fluctuate due to the different resonant structures, and |R − r| is the separation between the MM-charge Q and the MS-VALBOND (state-dependent charge) q.
For most of the structures along an MD trajectory the charge q k on a particular atom k will be that of one specific state i with weight w i ≈ 1 and all other weights w j ≈ 0. However, if two or several states are close in energy (see e.g. Figure 5B) or in the neighborhood of a crossing between two (or more) states (see Figure 10B) several weights will differ from 0. This is similar in the multi state adiabatic reactive MD approach where the weights, however, depend on the energy gap between the lowest energy state for a given conformation and all higher excited states. 49 Because the weights depend on the geometry of the part which is treated with MS-VALBOND, the charge on atom k will also depend on the geometry which effectively amounts to a fluctuating charge where q i k is the charge on atom k in state i and the total charge of the solute must be constant  numerically according to where, X and x are the coordinates of the solvent and solute atoms, respectively.

Dynamics of ClF 3 in Gas Phase and in Solution
For the parametrization of the diagonal states of ClF 3 the core fragment ClF + 2 was considered. A Morse potential was fitted to reference ab initio data for the Cl-F bond obtained at the CCSD/aug-cc-pVDZ 68

236.597
In the second step the parametrization of the off-diagonal terms was carried out by considering six degrees of freedom in ClF 3 system: Cl-F B , Cl-F C , Cl-F D , F B -Cl-F C , F C -Cl-F D and F B -Cl-F C (-)F D dihedral angle. The off-diagonal force field parameters for ClF 3 are also provided in Table 1.
With the FF available, first the system was energy minimized, heated to 300 K and then equilibrated for 500 ps. Further 500 ps of N V E simulations were performed for the analysis.
The solvated system was generated by immersing the solute molecule in a cubic TIP3P water box 70 of length 30Å .
As a validation, the energy minimized geometry using this parametrization was compared with experimental data 71 and with the optimized structure obtained at the CCSD/aug-cc-pVDZ level of theory (see Table 2). In general, the geometrical parameters for bond lengths and bond angles obtained from MS-VALBOND agree well with both sets of reference data.
The Cl-F B,D bond lengths from MS-VALBOND deviate by -0.02Å and -0.06Å from the CCSD calculations and experiment, respectively. Similar observations are made for the Cl-F C bond which differs by -0.04Å and -0.08Å from the CCSD calculations and experiment.
For the F B -Cl-F D angle the difference between MS-VALBOND calculations and the reference data is ∼ 2 • (see Table 2) which increases to 6 • for the F B -Cl-F C angle. All three approaches yield a planar geometry of the molecule.    To these references data two Morse potentials for Pt-N and Pt-Cl bonds were fitted and parameters are given in the SI (see Table S1).
For the parametrization of the off-diagonal terms, the element V 12 (transition between states I and II) is considered as an example. In going from resonance structure I to II, one Pt-N bond breaks and one Pt-Cl bond forms. Therefore, the off-diagonal term is described  Table 3. Cl Figure 6: The six resonance structures of cisplatin. Resonance structures II and IV are chemically equivalent to structures V and VI. NBO analysis with Natural Resonance Theory 39 demonstrates that all six resonance structures exhibit significant population depending on the molecular geometry.
For structural validation, the minimized geometry from MS-VALBOND was compared with the corresponding crystal structure, 73 see Figure 7 and Table 4. It is found that bond lengths and angles from MS-VALBOND agree to within less than 2 % with the X-ray data. Specifically, MS-VALBOND yields bond distances Pt-N and Pt-Cl in cisplatin of 2.07Å and 2.36 A, respectively, compared with 2.05Å and 2.32Å from the X-ray crystal structure. 73    A further validation of the MS-VALBOND force field for cisplatin was done by computing the infrared spectrum without adjusting the parametrization. The spectrum was calculated from the dipole moment auto-correlation function C(t) =< µ(0)µ(t) > from a 1 ns MD simulation in gas phase. The total dipole moment µ(t) was recorded along the trajectory and correlated over 2 15 time origins separated by 1 fs to yield C(t). Then the spectral signal C(ω) was computed from the fast Fourier transform (FFT) with a Blackman filter. Finally, the infrared spectrum A(ω) was determined by Boltzmann weighting C(ω) according to Here ω, k B and T are the frequency, Boltzmann constant and temperature, respectively. Figure S2 reports the computed IR-spectra (black trace) compared with the experimentally determined 80 vibrational frequencies (no intensity information is available from experiments).
A direct comparison for the most important modes is provided in Table 5 and shows that they compare favourably. This is particularly true for the L-M-L bending frequencies for all three N-Pt-N, Cl-Pt-Cl and N-Pt-Cl angles which match very well with their experimental frequencies 80 (see Table 5). On the other hand the Pt-N and Pt-Cl stretch frequency deviate by about 50 cm −1 from the experiments. 80 Finally, the peak at ∼ 900 cm −1 in the computed spectrum corresponds to the NH 3 wagging vibration and also agrees closely with experiment (at 927 cm −1 ). Further improvement could be achieved by explicitly fitting for example the stretching force constants to adjust the Pt-N, Pt-Cl or N-H stretch frequencies which was, however, not pursued here. Finally, the quality of the MS-VALBOND force field was also established by comparing to electronic structure calculations for thermally populated structures. For this an MD simulation at 50 K was carried out from which 50 structures were randomly selected. Their energies were determined at the B3LYP level of theory with the 6-31G(d,p) basis set for all atoms except Pt for which LANL2DZ basis is used. These energies are then correlated with the between DFT and MS-VALBOND is 0.33 kcal/mol which is within chemical accuracy and further demonstrates the quality of the MS-VALBOND parametrization for this TM complex.    halomethane substitution reactions at the MP2 and CCSD level of theory in implicit water also yielded favourable agreement with experiment whereas DFT underestimated the barrier height. 66,86 As a comparison, the gas-phase barrier for the BrMe + Cl − → Br − + MeCl reaction from CCSD/aug-cc-pVDZ calculations is found to be 11.7 kcal/mol, in close agreement with previously computed 11.2 kcal/mol and experimental values 9.5 kcal/mol. 87,88 The free energy profile for the gas phase reaction is reported in Figure S3.
For the MS-VALBOND parametrization of the BrMe + Cl − → Br − + MeCl S N 2 reaction, the diagonal states were described by two Morse potentials and fitted to reference data from CCSD/aug-cc-pVDZ calculations along the C-Br and C-Cl bonds in BrMe or MeCl, respectively. All the Cl/Br-C-H angles were described by VALBOND with a sp 3 hybridization.
For the off-diagonal states, the parametrization was carried out by considering two degrees To study the BrMe + Cl − → Br − +MeCl reaction in aqueous solution, two reactants Cl − and MeBr were solvated in a cubic TIP3P water box of length 30Å. The system is energy minimized, followed by heating and an equilibration of 1 ns at 300 K. The energy barrier of 24.7 kcal/mol (see above) associated with the substitution reaction BrMe + Cl − → Br − + MeCl is too high for directly sampling it in unbiased simulations on time scales accessible to conventional MD simulations. In order to validate the present implementation for following chemical reactions without using a bias the reaction barrier is first artificially reduced to 2.3 kcal/mol by decreasing the dissociation energy of the C-Br bond. With this modification the reaction occurs on the 10 ps time scale (see Figure 10) in unbiased MD simulations. The total energy is conserved to within 0.2 kcal/mol (see Figure 10) which confirms that the present implementation is also suitable to follow chemical reactions in a meaningful fashion. Next, the S N 2 reaction for BrMe and Cl − in aqueous solution was investigated using the true potential energy surface with a barrier height of 24.6 kcal/mol. For this, umbrella sampling 89 simulations were performed. The distance difference δ = d C−Br − d C−Cl was used as the reaction coordinate. The system was prepared as described above, followed by a total of 50 umbrella simulations along the reaction coordinate with −2.5 ≤ δ ≤ 2.5. For each umbrella window 250 ps of equilibrium simulations (100 ps equilibrium simulation and 150 ps for data accumulation) were carried out. For all umbrella windows a force constant k = 50 kcal/mol was used.
The weighted histogram analysis method (WHAM) 90 was used to determine the resultant free energy profile for the BrMe + Cl − → Br − + MeCl reaction in water, see Figure 11.
The computed free energy barrier is 24.6 kcal/mol which agrees well with the experimental 65,85 value of 24.7 kcal/mol. Also, the barrier for the reverse reaction is 26.7 kcal/mol which agrees quite well with experiment that finds 27.8 kcal/mol. 65,66 Finally, the computed free energy difference ∆∆G = 2.1 kcal/mol between reactant and product is also consistent with experiment which found 3.1 kcal/mol (see Figure 11 red line). 65,66 Probably, further improvement of the force field would be possible by slightly tuning force field parameters such as the C-H bonds and H-C-H angles. Figure 11: The free energy profile for the BrMe + Cl − → Br − + MeCl reaction in water obtained from the umbrella sampling simulations. Experimental values shown in red. 65 Atom color code used: C (cyan), Br (pink), Cl (green) and H (white). For the free energy profile in gas phase see Figure S3.

Conclusion
The multi-state VALBOND (MS-VALBOND) approach discussed in the present work has been successfully implemented and applied to a range of situations including hypervalent molecules, transition metal complexes and chemical reactions. It is demonstrated that total energy is conserved and, whenever possible, comparison with experiment leads to good agreement for structures, vibrational frequencies or free energy barriers. Hence, quantitative atomistic simulations for hypervalent and metal-containing systems on extended time scales become possible with the present implementation which allows to make contact with experiments and provide additional molecular-level insight from the analysis of the simulations as is possible from MS-ARMD studies. 23,91 This is one of the main differences compared with methods such as ReaxFF which are more suitable for qualitative investigations.
The current implementation is also suitable for applications to larger systems, such as proteins because the reactive subsystem (e.g. a protein active site) can be represented with MS-VALBOND whereas the remainder of the simulation is treated with an empirical force field. Given the flexible implementation it is expected that MS-VALBOND can be applied to a range of challenging problems, including transition metals in biological systems or reactions involving organometallic systems.