# data[].attributes.body
TITLE:
Computational methods to dissect gene regulatory networks in cancer.

AUTHORS:
Archana S. Iyer, Hatice U. Osmanbeyoglu, Christina S. Leslie.
DATE :
21 APRIL 2017

METHOD:
-Network inference method.
-Network diffusion method.
-MutSIG and GISTIC method.
-Integrative mutation / CNV analysis method.
- Cross- cutting methods.

WHAT THE RESEARCH IS ABOUT?
1. The research is about diverse computational methodologies that have sought to interpret somatic alterations and genes expression data through models of gene regulatory networks.
2. The development of these approaches have been divided into four categories.
a- network inference based on gene expression data.
b- Interpretation of genetic alterations with respect to a prior network.
c- Network integration of cis regulatory information with tumor gene expression.
d- Cross-cutting algorithms that link changes in upstream signaling, through a prior network or proteomics data, to downstream transcriptional programs and thus decipher the regulatory impact of somatic alterations.

URL :
https://www.sciencedirect.com/science/article/pii/S2452310017300793

# data[].attributes.body
TITLE:
Mg co-ordination with potential carcinogenic molecule acrylamide: Spectroscopic, computational and cytotoxicity studies

DATE AVAILABLE ONLINE: 27 October 2017
AUTHORS: Kamran Dastafkana, Armin Kianib, Ali Obeydavia, Mohammad Rahimic

COMPUTATIONAL METHOD:
- DFT method
- TD-DFT method
- Polarisable continuum model (PCM)

WHAT IS RESEARCH ABOUT?
1.  The research is about synthesis Acr-Mg complex and characterized it and identify the binding site of Acr and Mg.In order to identify the binding sites of acr with Mg, Mg complex of acr has been synthesized using magnesium nitrate and acr as the starting material.
2. three sample which are pure Acr,  magnesium nitrate and Acr-Mg is used to detect and measure the level of toxicity
3. In-vitro cytotoxicity studies revealed significant decrease in the toxic level of acr-Mg complex as compared to pure acr.
Cyto-toxicity studies of acr, magnesium nitrate and acr-magnesium (acr-Mg) complex were also carried out. teh Acr-Mg complex is less cyto-toxic than Acr and Magnesium nitrate when they are at same concentration
4.the toxicity studies have been also carried out by MTT assay against MDA-MB-231 human breast cancer cell line. the acr-Mg complex were observed with decreace in toxicity compared to pure acr.

url :
https://www.sciencedirect.com/science/article/pii/S0022286017314539

# data[].attributes.body
TITLE :

Structural and computational study on inhibitory compounds for endonuclease activity of influenza virus polymerase

AUTHORS :  Satoshi Fudo, Norio Yamamoto, Michiyoshi Nukaga, Takato Odagiri, Masato Tashiro, Saburo Neya Tyuji Hoshino
DATE AVAILABLE ONLINE :
29 July 2015
METHOD :
--Molecular dynamics simulation
--ONIOM method implemented in Gaussian09
--UFF
--A general AMBER force field (GAFF)
WHAT THE RESEARCH IS ABOUT:
1. The binding mode of one of the inhibitory compounds to PAN was investigated in detail by means of X-ray crystal structure analysis and molecular dynamics (MD) simulation.
It was observed in the crystal structure that three molecules of the same kind of inhibitor were bound to one PAN.
2. The last binding site was not observed in previous crystallographic studies. Hence, the stability of inhibitor binding was examined by performing 100-ns MD simulation.
3. During the MD simulation, the three inhibitor molecules fluctuated at the respective binding sites at different amplitudes, while all of the molecules maintained interactions with the protein.
4.  The changes in binding mode are consistent with the results of calculation for binding free energy, in which molecule 1C had the lowest value, and with the fact that electron density was most clearly observed for molecule 1C in the crystal structure.
URL : https://www.sciencedirect.com/science/article/pii/S0968089615006306

# data[].attributes.body
TITLE :

Size dependence of graphene chemistry: A computational study on CO desorption reaction

AUTHORS :
Buu Q. Pham, Vu H. Nguyen, Thanh N. Truong
DATE AVAILABLE ONLINE :
12 January 2016

METHOD :
Density Functional Theory
WHAT THE RESEARCH IS ABOUT :
1. The size-dependent chemistry of CO desorption from graphene ribbons was systematically studied using Density Functional Theory.
2. Two aspects of the size dependence were examined, namely the energy gap between the electronic low spin ground state and its high spin excited state and the fluctuation of the energy profiles of CO desorption reactions in the ground state.
3. To study the  effects of physical size on the relative barriers and stabilities of reaction intermediates along the CO desorption channel could be attributed to the disruption of π-conjugation of the graphene ribbon surface.
URL :
https://ezproxy.upm.edu.my:2060/science/article/pii/S0008622316300288

# data[].attributes.body
TITLE:

Experimental and computational study of membrane affinity for selected energetic compounds.

AUTHORS:
Golius A, Gorb L, Scott A.M, Hill F.C, Manoj Shukla, Goins A.B, Johnson D.R, Leszczynski J.
DATE AVAILABLE ONLINE:
25 January 2016
METHOD:
Quantitative Structure – Activity (Property) Relationships (QSAR/QSPR)
WHAT THE RESEARCH  IS ABOUT:
1. The affinity of various energetic compounds for a biological membrane was investigated using experimental and computational techniques. We measured octanol–water (log(Kow)) and liposome–water (log(Klipw)) partition coefficients for the following chemicals: trinitrotoluene (TNT), 2,4-dinitrotoluene (2,4-DNT), 2,6-dinitrotoluene (2,6-DNT), 2,4-dinitroanisole (DNAN), 2methoxy-5-nitrophenol (2M5NP), 2,4,6-trinitrobenzene (TNB), and 2,4-dinitrophenol (2,4-DNP). In order to determine log(Klipw) experimentally, we used artificial solid supported lipid liposomes produced under trade mark TRANSIL. Log(Kow) value were predicted with several program packages including the COSMOthermX software. Log(Klipw) were estimated with COSMOmic as implemented in the COSMOthermX program package.
2. To characterize the ability of a substance to bioaccumulate is to consider the partition coefficient between water and a lipid pool i.e. biological membrane. In its simplest form, bioconcentration has been characterized as the ratio of the equilibrium concentrations of a chemical (solute) that partitions between n-octanol (representing the lipid fraction) and water.
3. To improve the ease and repeatability of bioconcentration measurements, recently test membranes comprised of a class of artificial membrane vesicles made up from phospholipids (liposomes) were synthesized and made commercially available.
URL:
https://www.sciencedirect.com/science/article/pii/S0045653516300108

# data[].attributes.body
TITLE:

Computational chemistry study of solid and aqueous solution interface

AUTHOR:  Ugur Mart , Changho Jung , Michihisa Koyama , Momoji Kuboa , Akira Miyamoto
DATE AVAILABLE ONLINE: 7 January 2005
METHOD: Density functional theory and molecular dynamic
WHAT THE RESEARCH IS ABOUT:
1) Investigate and compare the ionic properties of monovalent metal ions on the hydration behavior and adsorption mechanisms at the solid–aqueous solution interface.
2) Understand how water molecules behave on an oxide surface in an electrolytic solution, and how the interatomic interactions compete to arrange the molecules.
3) An MgO cluster model, which includes nine Mg and nine O atoms, was used in order to investigate the interaction of hydrated cations with oxide surface.
4) Hydration of monovalent metal ions (Li+ , Na+ , K+ and Cs+ ) was carried out by aligning five water molecules initially around the cation to analyze how water molecules interact with single cation depending on its structure-making or structure-breaking behavior.
URL: https://ac.els-cdn.com/S0169433204018082/1-s2.0-S0169433204018082-main.pdf?_tid=f96c0d74-14b5-11e8-8d64-00000aacb35f&acdnat=1518963424_5a76ae0650e8fb57d8e2a1df2c7b36a7

# data[].attributes.body
Title :  New Optimization Method for Conformational Energy Calculations on Polypeptides: Conformational Space Annealing
Author :  Jooyoung L., Harold A.S., & Rackovsky S.
Date  : 8 January 1997
Method : Optimization method
What the research is about :
The novel aspect of the proposed CSA method is to keep track of many ‘‘distinct’’ local minima that are designed to be distributed as far apart as possible in the conformational space.  The CSA method maintains diversity in the bank by setting the value of D large cut enough at the beginning of the algorithm and slowly reducing it, hence conformational space annealing, and by selecting a new seed conformation significantly far away in the conformational space from the previous one.
URL :  http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.385.2677&rep=rep1&type=pdf

# data[].attributes.body
TITLE: A computational chemistry perspective on the current status and future direction of hepatitis B antiviral drug discovery
AUTHOR:  Morgnanesi D, Heinrichs EJ, Mele AR, Wilkinson S, Zhou S, Kulp JL 3rd.

DATE AVAILABLE ONLINE: 23 Oct 2015
WHAT THE RESEARCH ABOUT:
1) Three computational techniques most commonly used in HBV research are homology modeling, molecular docking, and molecular dynamics.
2) Homology modeling used to predict protein structure and to construct conformers of the viral polymerase and the HBV X protein.
3) Molecular docking is used to predict the most likely orientation of a ligand when it is bound to a protein, as well as determining an energy score of the docked conformation.
4) Molecular dynamics is a simulation that analyzes biomolecule motions and determines conformation and stability patterns.
5) All of these modeling techniques have aided in the understanding of resistance mutations on HBV non-nucleos(t)ide reverse-transcriptase inhibitor binding

METHODS: homology modeling, molecular docking, and molecular dynamics.
URL: https://www.ncbi.nlm.nih.gov/pubmed/26477294

# data[].attributes.body
MOEPGA: A novel method to detect protein complexes in yeast protein–protein interaction networks based on MultiObjective Evolutionary Programming Genetic Algorithm
Author: Buwen Caoa, Jiawei Luoa, Cheng Lianga, Shulin Wanga, Dan Song
Year: 2015
Date available online: 7 July 2015

Method: MultiObjective Evolutionary Programming Genetic Algorithm (MOEPGA)

What the research is about:

A novel MultiObjective Evolutionary Programming Genetic Algorithm (MOEPGA) which integrates multiple network topological features to detect biologically meaningful protein complexes. It is first systematically analyzes the multiobjective problem in terms of identifying protein complexes from PPI networks, and then constructs the objective function of the iterative algorithm based on three common topological properties of protein complexes from the benchmark dataset, finally describe algorithm, which mainly consists of three steps, population initialization, subgraph mutation and subgraph selection operation. To show the utility of this method, MOEPGA is compared with several state-of-the-art algorithms on two yeast PPI datasets. The experiment results demonstrate that the proposed method can not only find more protein complexes but also achieve higher accuracy in terms of fscore. Moreover, it cover a certain number of proteins in the input PPI network in terms of the normalized clustering score and  can serve as a powerful framework to detect protein complexes in yeast PPI networks, thereby facilitating the identification of the underlying biological functions
Link : https://ac.els-cdn.com/S1476927115300517/1-s2.0-S1476927115300517-main.pdf?_tid=b6b6a5d2-132c-11e8-9980-00000aacb35e&acdnat=1518794520_f41bb46d6e6de7430212894ec3d88573

# data[].attributes.body
Title: Mineral surface chemistry control for origin of prebiotic peptides
Author:Valentina Erastova,Matteo T.Degiacomi,Donald G.Fraser and H.Chris Greenwell
Date:11 December 2017
What the research about :
They carry out a large scale computational modeling study of interactions between amino acids and LDHs under reducing early Earth conditions. Their LDH layers have the composition of [Mg₃Al(OH)₈]⁺, while interacting with amino acids and peptides are deprotonated at pH 9.5. For this study, they chose a variety of amino acids (alanine, aspartate, leucine, lysine, histidine, and tyrosine), their mixtures, di- and hexa- peptides, and a randomly created 24-amino-acid-long peptides from the amino acid distributions mimicking naturally occurring ones. Besides,they explored the conceptual challenges associated with the formation of proto -biopolymers on mineral surfaces using the example of oligopeptides and amino acids at LDH surfaces.Some of the results obtained were intercalation of amino acids affects LDH layer dynamics, amino acids and peptides adsorb on LDHs via their C-termini,arrangement of adsorbed species is templated by LDH structure, LDHs promote amino acid polymerization ,peptide bond formation on LDHs is energetically favorable, and adsorbed amino acids and peptides diffuse on LDH surface.
Computational Methods:
1)Molecular dynamics simulations were performed with GROMACS 4.6.7
2)DBSCAN clustering algorithm
3)The radial distribution function (RDF) and the C-terminal atoms of the amino acids using LDH aluminum as reference.
4) The snapshots were produced with VMD 1.9.1 and graphs were produced with Matplotlib
5)The quantum calculations were performed with CASTEP
6)Developed a kinetic model to describe the process of peptide formation assisted by the LDH interlayer
7)Molecular models set-up

URL:
https://www.nature.com/articles/s41467-017-02248-y#Fig1

# data[].attributes.body
Computational methods to design cyclic peptides

By: Sean M McHugh, Julia R Rogers, Sarah A Solomon, Hongtao Yu, Yu-Shan Lin

1 September 2016

About: Cyclic peptides (CPs) are promising modulators of protein–protein interactions (PPIs), but their application remains challenging. It is currently difficult to predict the structures and bioavailability of CPs. The ability to design CPs using computer modeling would greatly facilitate the development of CPs as potent PPI modulators for fundamental studies and as potential therapeutics. Herein, computational methods are used to generate CP libraries for virtual screening, as well as current efforts to accurately predict the conformations adopted by CPs. These advances are making it possible to envision robust computational design of active CPs. However, unique properties of CPs pose significant challenges associated with sampling CP conformational space and accurately describing CP energetics. These major obstacles to structure prediction likely must be solved before robust design of active CPs can be reliably achieved.
Computational methods:
1)  Peplook, PEPstrMOD, PEP-FOLD, and I-TASSER ( determine conformations)

2) Molecular mechanics (to characterize the structural ensembles of CPs)

3) Macro Model (Macro Model, Schrodinger), MOE (MOE, Chemical Computing Group)

url: https://doi.org/10.1016/j.cbpa.2016.08.004

# data[].attributes.body
Title : In silico toxicology: computational methods for the prediction of chemical toxicity
Year : 6 January 2016

Author : Arwa B. Raeis, Vladimir B. Bajic
Method :  in silico toxicologymethods
What the research about :
1)Determining the toxicity of chemicals to identify their harmful effects on humans, animals, plants, or the environment
2)In silico toxicology aims to complement existing toxicity tests to predict toxicity, prioritize chemicals, guide toxicity tests, and minimize late-stage failures in drugs design
3)The Tox-21c stresses replacing animal testing with human-relevant testing methods, either in vitro or in silico
url :  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4785608/

# data[].attributes.body
Title: Prediction of quartz crystal microbalance gas sensor responses using a computational chemistry method
Author: K.Nakamura, T.Nakamoto and T.Moriizumi
Date: 14 December 1999
What the research is about:
The research is for predicting the quartz crystal microbalance (QCM) gas sensor responses using computational chemistry method and the grand canonical Monte Carlo (GCMC) method. They also tried to calculate the amounts of sorption of gas molecules such as alcohols, aromatics, ketones, esters, alkanes and perfumes into the sensing films. Next, the research is to compare the predicted partition coefficients of alcohols, aromatics, ketones, esters, alkanes and perfumes for typical GC films with the experimental ones.

Computational Methods:
1) Using Cerius2 (BIOSYM/MSI, Cerius2) for the simulation
2) The computer used for the calculation was Indigo2, Silicon Graphics
3) The GCMC (grand canonical Monte Carlo) method was used for the sorption simulation
URL: https://doi.org/10.1016/S0925-4005(99)00071-4

# data[].attributes.body
184743

A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP)

Takeshi Yanai, David P Tew, Nicholas C Handy
21 July 2004

-A new hybrid exchange–correlation functional named CAM-B3LYP is proposed. It combines the hybrid qualities of B3LYP and the long-range correction presented by Tawada et al. [J. Chem. Phys., in press]. We demonstrate that CAM-B3LYP yields atomization energies of similar quality to those from B3LYP, while also performing well for charge transfer excitations in a dipeptide model, which B3LYP underestimates enormously. The CAM-B3LYP functional comprises of 0.19 Hartree–Fock (HF) plus 0.81 Becke 1988 (B88) exchange interaction at short-range, and 0.65 HF plus 0.35 B88 at long-range. The intermediate region is smoothly described through the standard error function with parameter 0.33

-Coulomb-attenuating method

-URL :
https://www.sciencedirect.com/science/article/pii/S0009261404008620

# data[].attributes.body
Title: Investigation of anticancer properties of caffeinated complexes via computational chemistry methods

Author: Koray Sayin, Ayhan Ungordu

Date available online: 6th Dec 2017

Method used:
1) Gauge-including-atomic-orbital (GIAO) method in vacuo. [B3LYP]
2) Density Functional Theory (DFT) used to get Molecular Electrostatic Potential (MEP) map, MEP contour and Frontier Molecular Orbitals (FMOs)
3) Molecular Docking method in water phase.

What the research is about:
The research is to investigate the structural, electronical and biological properties of theobromine and their tautomers and Ruthenium and Osmium complexes. This is to study their anti-cancer activity. This research is done to find out about other metals than Pt to be used as drug to treat cancer.

Photos: https://drive.google.com/drive/folders/13H88f3lZt-NJmTWARkraPUyWDgaU8w6u?usp=sharing

# data[].attributes.attachment_link.title
Investigation of anticancer properties of caffeinated complexes via computational chemistry methods

# data[].attributes.body
Ab initio study of ionic solutions by a polarizable continuum dielectric model

Maurizio Cossia, Vincenzo Baronea, Benedetta Mennuccib, Jacopo Tomasib

10 April 1998

What the research is about :
- new implementation of integral equation formalism version of the PCM solvation method
- this implementation number of techniques for analysis of energies and electronic properties of closed and open shell molecule at the Hartree-Fock and density functional levels
- HF calculations was implemented on closed shell systems and the present formulation is more effective then extended to DF calculations for open shell molecule

Computational method
- the solute molecule is embedded in a cavity surrounded by a polarizable medium, with the dielectric constant of the considered solvent. The cavity has molecular shape: it is formed by fused spheres, centered on solute atoms or atomic groups, according to the GEPOL algorithm
- IEF-PCM exploited in theory of integral equations
- IEF differs from the other PCM-like approaches only for Gel, in particular, in ionic solutions the ionic atmosphere affects the electrostatic solute–solvent interactions only

URL
https://www.sciencedirect.com/science/article/pii/S0009261498001067?via%3Dihub

# data[].attributes.body
Reliable Charge Assessment On Encapsulated Fragment For Endohedral Systems

A. J. Stasyuk, M. Solà & A. A. Voityuk
(2018)

- To determine charge distribution in endohedral complexes.

- Present a new scheme for assessment of charge distribution on endohedral complexes
- All calculations were performed with the PBE0/Def2-TZVP (with ECP-28 and ECP-60 for Ag and Au atoms³⁹ correspondently) method⁴⁰,⁴¹,⁴² using Gaussian 09 program⁴³.

- Charge density analysis performed with Mulliken, Löwdin, Hirshfeld, CM5, QTAIM, and NPA schemes for endohedral X@B39 complex and the predicted charge of X (Qx predicted)a obtained within the method proposed in this work. Units are electrons.

# data[].attributes.attachment_link.readable_text
Altmetric: 2

More detail

Article | Open

Scientific Reports 8, Article number: 2882 (2018)

doi:10.1038/s41598-018-21240-0

Download Citation

Received:

Accepted:

Published online:

Abstract

A simple scheme to determine charge distribution in endohedral complexes is suggested. It is based on comparison of inner-shell atomic orbital energies of the encapsulated species to the corresponding energies in reference systems with unambiguously defined charges on X. This robust approach is applied to endohedral borospherenes X@B ₃₉, for which the conventional schemes provide in some cases quite different results. Efficiency of proposed scheme also has been proven for typical fullerene based Sc ₃N@C ₈₀ endohedral complex.

Introduction

Atomic charge is one of the most widely used concepts in chemistry. But, it can unambiguously be defined only in some trivial cases. For instance, in homogeneous diatomic molecules with the total charge Q, the charge on each atom is equal to Q/2. In general, however, some assumptions should be made in order to quantify the charge distribution in a molecular system. Many different schemes to define atomic charges have been developed. Most popular of them are schemes based on the distribution of the electron density obtained from quantum-mechanical calculation and its subsequent transformation into atomic charges -Mulliken1, Lowdin2, Hirshfeld3, CM54, Weinhold5,6 (natural population analysis, NPA), and Bader7,8 (quantum theory of atoms in molecules, QTAIM) population analysis.

Discovery of buckminsterfullerene C ₆₀ 9, and its La encapsulated derivative10,11,12 made a new milestone in chemistry. Fullerenes containing confined atoms or small molecules are called endohedral fullerenes (EFs). Most often the incarcerated species is a metal atom or metallic cluster, and such systems are called endohedral metallofullerenes (EMFs)13,14,15,16,17,18,19,20,21,22. Charge transfer between the encapsulated species and the cage in EMFs is important because it determines the most suitable cage among the many possible isomers. For instance, in Sc ₂C ₂@C ₆₈ with a formal charge transfer of 4 electrons, the cage corresponds to the 6073 isomer23, whereas the cage that encapsulates Sc ₃N@C ₆₈ with a formal charge of 6 electrons is the 614024. Although many EMFs have been intensively studied they continue to attract much attention of researchers from various fields of natural sciences. Recent discovery by Zhai et al. of all-boron fullerenes B ₄₀− and B ₄₀ ⁰ (referred in literature as borospherenes)25 has given a new impetus to the chemistry of endohedral compounds and stimulated a large number of experimental and theoretical studies. A few months later Bai et al. studied computationally the viability of endohedral metalloborospherenes M@B ₄₀ (M = Ca, Sr). DFT calculations revealed great stability due to the almost perfect match in size of encapsulated fragment and cage. It was found that all these metalloborospherenes are formally charge-transfer M ²⁺B ₄₀ 2− complexes26. Other endohedral complexes of borospherenes based on B ₃₆− 27, B ₃₈ 2− 28, and B ₃₉− 29, cages demonstrate similar behaviour. A number of works predicting the existence of endohedral borospherenes have been published very recently30,31,32,33,34. Making an analogy with carbon fullerenes we can expect that endohedral borospherenes will be synthetized at an early date.

Results and Discussions

Focusing our attention on endohedral Cl@B ₃₉ borospherene complex based on B ₃₉ cage subunit, we have encountered the problem of the reliability of charge determination on central fragment. The population analysis carried out within various schemes demonstrates significantly different distribution of charges. However, a reliable definition of the charge localized on the encapsulated fragment as well as its changes is important for the understanding of chemical reactivity, and, in particular, for photo-induced charge transfer (CT) process in endohedral compounds35,36,37,38.

In this work, we present a new scheme for assessment of charge distribution on endohedral complexes. Proposed method is based on straightforward physical model and could be used as a convenient tool for the culling of inadequate charge schemes. All calculations were performed with the PBE0/Def2-TZVP (with ECP-28 and ECP-60 for Ag and Au atoms39 correspondently) method40,41,42 using Gaussian 09 program43.

Numerous proofs that metal-encapsulating borospherenes could demonstrate remarkable charge transfer properties have been recently reported. Among them, the work published by Chen et al. 29 demonstrates at first-principles level the viability of the axially chiral metalloborospherenes Ca@B ₃₉⁺. At the same time there are no examples of halogen encapsulated borospherenes, despite the fact that such complex is a unique system of halogen encapsulate into superhalogen unit. The uniqueness of the Cl@B ₃₉ complex is that both fragments exhibit high electron affinity (Cl - 3.612 eV44, B ₃₉ - 3.845 eV45). Thus it seems extremely difficult to predict charge distribution in such system. With a certain degree of confidence it can be assumed that charge separation should be observed, however it is completely unclear which of the fragments will be negatively charged and which positively.

Charge distribution analysis for B₃₉ boroshperene complexes

We focused our attention on complexes such as Cl@B ₃₉ (Fig. 1a), in particular on the description of charge states in such systems using various approaches. We found that charges on internal fragment, as well as on the cage, depend dramatically on the population scheme used.

Figure 1

Graphical representation of endohedral boroshperene complex Cl@B ₃₉ ( a) and (B ₃₉ + Cl) vdW complex ( b).

In contrast, for van der Waals complex (Fig. 1b) composed of B ₃₉ and Cl separated by 3 Å, (B ₃₉ + Cl) vdW, the predicted charges demonstrate a significant similarity (Table 1).

Table 1: Charge density analysis (in electrons) performed with Mulliken, Löwdin, Hirshfeld, CM5, QTAIM, and Natural Population Analysis (NPA) schemes for endohedral Cl@B₃₉ and (B₃₉ + Cl) van der Waals (vdW) complex obtained at PBE0/Def2-TZVP level of theory. Total charge is in all cases zero.

As can be seen from Table 1, for Cl@B ₃₉ complex, charge on the Cl fragment can vary from −0.62 e (QTAIM) to 0.76 e (Löwdin) depending on the used scheme. At the same time, charge on the cage ranges from −0.76 e (Löwdin) to 0.66 e (QTAIM). Presuming that behaviour of these fragments in endohedral complex does not change significantly, the reliability of the charges obtained based on Löwdin population analysis ( Q Cl = 0.76), as well as on Hirshfeld ( Q Cl = 0.00) and CM5 ( Q Cl = 0.14) schemes should be called into question. For the van de Waals complex, the obtained charge values on the Cl atom are much closer to each other and vary from −0.24 e (Löwdin) to −0.58 e (QTAIM).

The main problem in this case is the selection of the most reliable method/scheme of charge assessment, which is almost impossible without resorting to chemical intuition.

Model description

It is well known that the ionization potential (binding energy) of the inner-shell electrons of an atom in a molecular system depends on the charge state and chemical surroundings of the atom46. A linear relationship of the binding energy (E b) measured by X-ray photoelectron spectroscopy (XPS) and an atomic charge q computed with quantum mechanical methods (Eq. ( 1)) was established for different classes of chemical compounds

The XPS shift ΔE b is determined with respect to the binding energy E b0 found for a reference. Further studies showed, however, that the electrostatic potential φ created by chemical surroundings has to be accounted for47

There are inherent physical limitations that prevent accurate determination of the atomic charges. In particular, the parameters q and φ in Eq. ( 2) refer to the ground electronic state of the system, while the measured XPS shift, ΔE b, reflects both the ground and the final excited state (with a hole in the inner shell). To separate the contributions Aq and φ, some additional assumptions should be used. In spite of these and some other restrictions, XPS has been widely employed for chemical analysis48.

Let's consider a model system consisting of a multi-electron fragment X with partial charge Q ref incorporated into a cage constructed of N atoms with partial charge Q ext (Fig. 2).

Figure 2

Physical model used in this work.

Similar to Eq. ( 2) we assume that the orbital energy of core electrons of an atom of fragment X (X consists of either a single atom or several atoms) can be approximately presented as a simple function of the charge Q X on fragment X, corrected by the external electrostatic potential φ ex on atom due to the cage

The non-trivial point here is how to calibrate ε k c o r r ( , ϕ e x) so that the charge Q X in the reference systems is uniquely defined.

We propose the following step-by-step scheme to derive the charge on X in X@Cage species:

1. 1.
One constructs reference systems where charge on X is unambiguously defined. The reference systems for such a situation include the fragment X carrying the total charge Q ref (e.g. Q ref = −1, 0, + 1 etc) embedded into a sphere having the charge Q ext. The total charge Q tot of the reference system, Q tot = Q ref + Q ext, is equal to the charge of X@Cage. In most cases X@Cage is neutral (Q tot = 0) and thus for any reference system Q ext = −Q ref. If X@Cage is charged, Q ext = Q tot − Q X. The radius of the sphere R is determined by the cage geometry. For a cage comprising N atoms, R can be defined as

1 R= 1 N ∑ i=1 N 1 r i

where r i is the distance from atom i in the cage to the center of the cage determined by the non- weighted mean of the coordinates of the atoms in the cage. In this approach:

2. 2.
Quantum mechanical calculation with molecular charge Q ref (e.g. Q ref = −1, 0, + 1, etc.) are carried out. The computed orbital energy of inner shell K, ( Q r e f), together with the electrostatic potential created by the sphere provides :

ε K c o r r ()= ε K ()+ R

The found values are used to obtain a quadratic interpolation of the dependence. Note that the encapsulated species X comprise usually only several atoms and its calculation is very fast.

3. 3.
Quantum mechanical calculations of X@Cage provide the value in the "real" system. Using the interpolated function derived in the previous step we could directly obtain the charge Q X on the encapsulated fragment.

These three steps represent the simplest scheme to estimate charge separation between encapsulated species X and the cage in endohedral fullerenes and in similar structures. As shown below, a more elaborated treatment of the electrostatic embedding in the reference systems does not provide significantly different results. Also, as we will see, the proposed scheme is quite robust to choosing the inner shell K in X when constructing the interpolated function ε k c o r r ( Q r e f).

The classical point charge model was used to evaluate the electrostatic correction. This model implies the replacement of cage atoms by array of point charges placed at the nuclei of the constituted atoms. This approach has proven itself well in various fields of chemical science and computational biology49,50.

To verify the applicability of the proposed model to the charge assessment of the central fragment in endohedral complexes, the borospherene species Cl@B ₃₉ was considered in detail (Fig. 3).

Figure 3

Transformation of Cl@B ₃₉ complex within the point charge model framework necessary for calibration curve construction. Left part of figure corresponds to Cl@B ₃₉ complex, right part represents central fragment (Cl) surrounded by compensation charges (black spheres) allocated at former boron nuclei places.

Primary, it is necessary to perform a series of calculations for Cl fragment in various charge states. For this we have selected 5 states - Cl 2−, Cl 1−, Cl ⁰, Cl ¹⁺, and Cl ²⁺ in their ground states. Despite the fact that Cl ²⁺ as well as other halogen cations goes beyond conventional chemical principles and its existence is apparently impossible under ordinary conditions, we will take it into account as a model point to have a wider range of calibration. Including a large set of charge states into calibration makes a proposed scheme more robust and universal due to the fact that even in difficult cases a result can be found within the interpolation grid. Of course, in case of Cl fragment, data about its possible charge states is well known, but in case of more sophisticated encapsulated unit a necessary information can simply not be available. Moreover, it was found that inclusion of 1s orbital energy of Cl ²⁺ unit does not lead to physically meaningless result, therefore we recommend to use greater number of charge states for constructing a calibration curve.

To take into account correction to electrostatic potential, boron atoms forming cage have to be replaced with point charges (Fig. 3) in such way that total compensation charge has to be equal to the charge on Cl fragment in current charge state. For example, for Cl ¹⁺ cation compensation charge have to be equal to −1, while for Cl 1− fragment the compensation charge have to be equal to +1. It seems obvious that there are more than one way of the compensation charge distribution over cage nodes when the neutrality condition is satisfied.

To investigate the effect of compensation charge distribution on 1s orbital energy we performed calculations for 4 different models (Table S1, ESI). To construct first of them, even charge distribution over the cage were used, while for Cl ⁰ charge state no compensation point charges have been applied (CCS1, Table S1). In next three schemes, the charges on cage were generated randomly to compensate the charge −1 and then be linearly scaled to the charge state −2. In the case of Cl ¹⁺ and Cl ²⁺, compensation charges on the cage were obtained from previously generated only with opposite sign (CCS2, CCS3, and CCS4 in Table S 1). Last scheme is characterized by random generation of charges on all cage nodes for each state of Cl fragment. Only equality of the charge on Cl fragment and sum of compensation charges was monitored (CCS5, Table S1). Calibration curves and equations of approximate functions are shown in Fig. S 1. Detailed analysis of given charge distribution schemes revealed that electrostatic correction for 1 s orbital energy value depends only very weakly on particular compensation charge scheme. For the studied charge distributions, the range of the correction values is about one order of magnitude less than the correction values themselves (Table S2). At the same time, central fragment charges predicted on the basis of the mentioned charge compensation schemes are in the range from −0.44 e to −0.51 e. Thus, taking into account small differences in charge prediction caused by particular charge compensation method, we will use the even charge distribution as a compensation method.

At the next stage of our work, we studied the effect of the charge states taken into consideration on the predicted charge values of the Cl@B ₃₉ complex. It was found that the charge predicted within the proposed method depends weakly on charge states of Cl fragments that were taken into account in calibration. Four different calibration curves were built using various Cl charge states - 1  st: [Cl 2−...Cl ²⁺] with 5 reference points ( Q Cl = −0.475), 2 nd: [Cl 2−...Cl ¹⁺] with 4 points ( Q Cl = −0.473), 3 rd: [Cl 1−...Cl ²⁺] with 4 points ( Q Cl = −0.470), and 4 th: [Cl 1−...Cl ¹⁺] with 3 points ( Q Cl = −0.470) (Fig. 4).

Figure 4

Calibration curves constructed using various number of reference points (top left - [Cl 2−...Cl ²⁺] - 5 points, top right - [Cl 2−...Cl ¹⁺] - 4 points, bottom left - [Cl 1−...Cl ²⁺] - 4 points and bottom right - [Cl 1−...Cl ¹⁺] - 3 points) for charge prediction on central Cl fragment in Cl@B ₃₉ complex (Cl 1s orbital energy in complex of interest = −101.628 a.u.).

As can be seen in Fig. 4, minor differences (of less than 0.01 e) in the predicted charge on the Cl fragment obtained from different calibration curves are found. This error is one order of magnitude smaller than the error caused by the use of different charge compensation schemes.

Thus, we can conclude that for studied Cl@B ₃₉ complex proposed method has shown its best side. However, we were wondering if proposed method can equally work well with other K-type orbitals. To check this, Cu@B ₃₉ complex has been considered. Cu atom has 3 inner s-orbitals. Calibrations curves construct for 1s, 2s, and 3s orbitals based on five Cu charge states [Cu 2−...Cu ²⁺] reference points predict charge on central fragment 0.416 e, 0.391 e, and 0.398 e correspondently (Fig. S 2, ESI). Thus, one can assume that any inner orbital, as far as it not involved in any orbital interaction with the cage orbitals, could be used for construction of the calibration. The latest suggests that the proposed scheme can be extended to "heavy" elements, for which pseudopotentials are used.

Model verification

Finally, to verify the proposed approach several systems with general formula X@B ₃₉, where X represents halogens (F, Cl, Br), small radicals (NO, CN, FO, CF) or metals of 1B group of periodic chart (Cu, Ag, Au) were investigated (Fig. S 3). The summary results are presented in Table 2.

Table 2: Charge density analysis performed with Mulliken, Löwdin, Hirshfeld, CM5, QTAIM, and NPA schemes for endohedral X@B₃₉ complex and the predicted charge of X (Q x predicted)a obtained within the method proposed in this work. Units are electrons.

As can be seen from Table 2, the proposed charge assessment approach provides the physically adequate results in all studied cases. For each cases, Mulliken and Löwdin charge assessment schemes provide the worst performance with correlations with the Q X predicted charges of R ² = 0.388 for Mulliken and R ² = 0.172 for Löwdin, respectively. Moreover, sometimes they provide physically meaningless results, such as in case of Au@B ₃₉ (Mulliken) or Cl@B ₃₉, Br@B ₃₉, FO@B ₃₉ and CN@B ₃₉ (Löwdin) complexes. Hirshfeld charges are quite well correlated with the predicted charges ( R ² = 0.909) but the charges are noticeably underestimated. CM5 method yields comparable to Hirshfeld correlation with the predicted results with R ² = 0.896. Equally to Hirshfeld scheme, significant charge underestimation for CM5 is observed. In some cases Hirshfeld and CM5 schemes does not reproduce the sign of charge on the central fragment (Fig. S 4, Table S3, ESI). Finally, the best correlations between calculated and predicted charges were found for NPA and QTAIM schemes (Fig. 5).

Figure 5

Correlations of the NPA and QTAIM charges with the results obtained by the proposed method. Units are electrons.

Thus, we can conclude that the predicted charge assessment method provide results comparable with well-known methods and can be used as a fast and convenient selection tool that avoids assignment of unphysical charges to encapsulated species.

To demonstrate that the proposed method can be applied not only to boron cages but also to a wide range of endohedral compounds, one of the most studied (both synthetically51,52,53 and theoretically54,55) Sc ₃N@I h-C ₈₀ cluster has been investigated. Detailed charge distribution analysis within Mulliken, Löwdin, Hirshfeld, CM5, QTAIM, and NPA schemes has been performed. For Sc ₃N@I h-C ₈₀ complex, calculated charge values on the Sc ₃ N fragment vary from −0.090 e obtained with Mulliken scheme to 0.851 e and 0.969 e for Löwdin and Hirshfeld ones, correspondently. CM5 scheme provides 2.062 e, while NPA and Bader analysis result in 3.426 e and 3.772 e (for more details see ESI). Scheme proposed in current paper predicts the charge on central fragment of 1.502 e based on 1s N orbitals calibration. Thus, with this latter example we show that encapsulated fragment charge estimation scheme described above can be universally used for any endohedral system.

The remarkable feature of this approach is that for fragment X consisting of more than one atom, calibration curve can be constructed using the 1s orbital energy of any of the atoms of the incarcerated molecule. For example, calibration curves constructed using O 1s orbital energies and N 1s orbital energies for NO@B ₃₉ complex predict 0.212 e and 0.253 e charge on NO species, respectively (in detail, the procedure for constructing the calibration curves and the charges prediction for denoted NO species and all others considered in this paper is presented in supporting information). An agreement between the results obtained from different calibration curves is a good indicator of reliability of the result. However, it should be emphasized that this charge assessment method does not allow determination of charge on the individual atoms of the encapsulated fragment X, but only provides the charge on the whole fragment. On the other hand, taking into account physical model underlying the method described the following limitations must be noted. First, the proposed model does not take into account orbital interactions. Thus, if orbitals of fragment X are involved in strong orbital interactions, a significant error could be observed. Second, because the compensation charge in our method were treated within the framework of point charge model, a small underestimation of the predicted charge could take a place. The reason for this is that charge associated with particular node of cell is not point charge in fact, but only distributed electron density located somewhat closer to the fragment X, which in turn, will certainly affect the accuracy of the electrostatic correction.

Conclusions

In conclusion, we have developed a new approach for charge assessment of the encapsulated fragment based on the comparison of the low-lying orbital energies of atoms in central fragment in complex of interest with reference systems. The proposed approach demonstrates excellent performance on endohedral borospherenes X@B ₃₉ with encapsulated metal atoms, halogens or small radicals. Moreover, the workability of proposed approach has been demonstrated on the typical fullerene based Sc ₃N@I h-C ₈₀ endohedral complex. In view of the physical and computational simplicity, the proposed method can be applied to very large systems. In cases when conventional schemes provide essentially different results, the proposed method could be used as a convenient and robust tool to exclude unreliable data from consideration. We are convinced that the proposed approach will find application in research of all kinds of endohedral complexes and related fields of chemical science.

Additional information

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. 7.
Bader, R. F. W. Atoms in Molecules: A Quantum Theory (Oxford, 1990).

2. 9.
Kroto, H. W., Heath, J. R., O'Brien, S. C., Curl, R. F. & Smalley, R. E. C ₆₀: Buckminsterfullerene. Nature 318, 162-163, https://doi.org/10.1038/318162a0 (1985).

3. 12.
Cox, D. M., Trevor, D. J., Reichmann, K. C. & Kaldor, A. C ₆₀La: a deflated soccer ball? J. Am. Chem. Soc. 108, 2457-2458, https://doi.org/10.1021/ja00269a060 (1986).

4. 20.

5. 22.
Popov, A. A. Endohedral Fullerenes: Electron Transfer and Spin (Cham, 2017).

6. 23.
Shi, Z.-Q., Wu, X., Wang, C.-R., Lu, X. & Shinohara, H. Isolation and Characterization of Sc ₂C ₂@C ₆₈: A Metal-Carbide Endofullerene with a Non-IPR Carbon Cage. Angew. Chem. Int. Ed. 45, 2107-2111, https://doi.org/10.1002/anie.200503705 (2006).

7. 24.

8. 27.
Tian, W.-J. et al. Saturn-like charge-transfer complexes Li ₄&B ₃₆, Li ₅&B ₃₆⁺, and Li ₆&B ₃₆ ²⁺: exohedral metalloborospherenes with a perfect cage-like B ₃₆ 4− core. Phys. Chem. Chem. Phys. 18, 9922-9926, https://doi.org/10.1039/C6CP01279E (2016).

9. 28.
Chen, Q. et al. Endohedral Ca@B 38: stabilization of a B ₃₈ 2− borospherene dianion by metal encapsulation. Phys. Chem. Chem. Phys. 18, 11610-11615, https://doi.org/10.1039/c5cp06169e (2016).

10. 29.
Chen, Q. et al. Endohedral C ₃Ca@B ₃₉⁺ and C ₂Ca@B ₃₉⁺: axially chiral metalloborospherenes based on B ₃₉. Phys. Chem. Chem. Phys. 17, 19690-19694, https://doi.org/10.1039/C5CP03178H (2015).

11. 31.
Li, S. X., Zhang, Z. P., Long, Z. W. & Qin, S. J. Structures, stabilities and spectral properties of borospherene B ₄₄− and metalloborospherenes MB ₄₄ 0/− (M = Li, Na, and K). Sci. Rep. 7, 40081, https://doi.org/10.1038/srep40081 (2017).

12. 32.
Xu, Q. H. et al. Computational investigation on MB n (M = Li-Cs, Be-Ba, Sc-La and Ti; n = 28 and 38). J. Mol. Model. 22, 184, https://doi.org/10.1007/s00894-016-3055-4 (2016).

13. 33.
Bin, L. et al. B ₄₀ fullerene as a highly sensitive molecular device for NH ₃ detection at low bias: a first-principles study. Nanotechnology 27, 075501, https://doi.org/10.1088/0957-4484/27/7/075501 (2016).

14. 35.
Kawashima, Y., Ohkubo, K. & Fukuzumi, S. Efficient Charge Separation in Li⁺@C ₆₀ Supramolecular Complexes with Electron Donors. Chem. - Asian J. 10, 44-54, https://doi.org/10.1002/asia.201403075 (2015).

15. 36.
Supur, M. et al. Graphene oxide- Li⁺@C ₆₀ donor-acceptor composites for photoenergy conversion. Phys. Chem. Chem. Phys. 17, 15732-15738, https://doi.org/10.1039/c5cp01403d (2015).

16. 43.
Gaussian 09 Rev. E.01. (full citation in ESI).

17. 44.

18. 46.
Slegbahn, K. et al. ESCA Atomic, Molecular and Solid State Structure Studied by Means of Electron Spectroscopy (Uppsala, 1967).

19. 47.
Slegbahn, K. et al. ESCA Applied to Free Molecules (Amsterdam/London, 1969).

20. 49.

21. 50.
Tsunekawa, S., Sahara, R., Kawazoe, Y. & Kasuya, A. Origin of the blue shift in ultraviolet absorption spectra of nanocrystalline CeO 2−x particles. Mater. Trans., JIM 41, 1104-1107 (2000).

22. 52.
Iezzi, E. B. et al. A Symmetric Derivative of the Trimetallic Nitride Endohedral Metallofullerene, Sc ₃N@C ₈₀. J. Am. Chem. Soc. 124, 524-525, https://doi.org/10.1021/ja0171005 (2002).

23. 55.
Aroua, S., Garcia-Borràs, M., Osuna, S. & Yamakoshi, Y. Essential Factors for Control of the Equilibrium in the Reversible Rearrangement of M ₃N@I h-C ₈₀ Fulleropyrrolidines: Exohedral Functional Groups versus Endohedral Metal Clusters. Chem. - Eur. J. 20, 14032-14039, https://doi.org/10.1002/chem.201403743 (2014).

Download references

Acknowledgements

We are grateful for financial support from the Spanish MINECO (CTQ2017-85341-P and CTQ2015-69363-P projects), the Catalan DIUE (2014SGR931, XRQTC, and ICREA Academia 2014 Award to M.S.), and the FEDER fund (UNGI10-4E-801). A. J. S. gratefully acknowledges The Interdisciplinary Centre for Mathematical and Molecular Modelling of the University of Warsaw (ICM) for computational facilities (grant no. G33-17).

Author information

Affiliations

1. Institut de Química Computacional and Departament de Química, Universitat de Girona, C/ Maria Aurèlia Capmany 69, 17003, Girona, Catalonia, Spain

A. J. Stasyuk

, M. Solà

& A. A. Voityuk

2. Institució Catalana de Recerca i Estudis Avancats, 08010, Barcelona, Spain

Authors

1. Search for A. J. Stasyuk in:

2. Search for M. Solà in:

3. Search for A. A. Voityuk in:

Contributions

A.J.S., M.S. and A.A.V. conceived and designed research, analyzed data, interpreted results and wrote the paper. A.J.S. performed the calculations.

Competing Interests

The authors declare no competing interests.

Corresponding authors

Correspondence to A. J. Stasyuk or M. Solà or A. A. Voityuk.

Electronic supplementary material

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

# data[].attributes.attachment_link.title
Reliable charge assessment on encapsulated fragment for endohedral systems

# data[].attributes.attachment_link.description
Altmetric: 2 More detail Article | Open Scientific Reports 8, Article number: 2882 (2018) doi:10.1038/s41598-018-21240-0 Download Citation Received: Accepted: Published online: A simple scheme to determine charge distribution in endohedral complexes is suggested. It is based on comparison of inner-shell atomic orbital energies of the encapsulated species to the corresponding energies in reference systems with unambiguously defined charges on X.

# data[].attributes.body
Applications of computational chemistry to the study of the antiradical activity of carotenoids: A review

Debora Luana Monego, Marcelo Barcellos da Rosa, Paulo Cicero do Nascimento
15 February 2017

URL:
https://www.sciencedirect.com/science/article/pii/S0308814616313139

- Most computational studies apply to the electron transfer mechanism. A significant number of these applications have been directed towards understanding the electron transfer mechanism, and a useful tool called the FEDAM (full-electron donor-acceptor map) was developed to better evaluate this mechanism.  Hydrogen abstraction may be favored according to the radical’s nature.

# data[].attributes.headline
Nurul Syazwani Azyan Zairul Akmal 182192

# data[].attributes.subject
Nurul Syazwani Azyan Zairul Akmal 182192

# data[].attributes.body
Computational insights for the hydride transfer and distinctive roles of key residues in cholesterol oxidase

Li-Juan Yu, Emily Golden, Nanhao Chen, Yuan Zhao, Alice Vrielink & Amir Karton (2017)

A number of experimental studies investigated the structure and catalytic mechanism of ChOx over the past two decades. The active site of ChOx is shown. With regard to the enzymatic mechanism, Sampson et al. identified three residues (Glu361, His447, and Asn485) that are significant for the enzymatic activity and Vrielink et al.⁴,⁵ revealed that Gly120 also plays an important role in the catalytic mechanism.
The method used is the QM/MM MD simulations. Development and applications of combined QM/MM methods for enzymes have been extensively reviewed in the past. In order to obtain the barrier height of the hydride transfer, QM/MM MD simulations have been carried out for this enzyme, considering the influence from the surrounding residues to the reaction site. Then, high-level double-hybrid DFT simulations. High-level DHDFT calculations using the B2GP-PLYP functional were performed in order to obtain accurate reaction energies and barrier heights of the hydride transfer from the C3 atom of the substrate to the N5 atom of the FAD cofactor.

# data[].attributes.attachment_link.readable_text
Altmetric: 7

More detail

Article | Open

Scientific Reports 7, Article number: 17265 (2017)

doi:10.1038/s41598-017-17503-x

Download Citation

Received:

Accepted:

Published online:

Abstract

Cholesterol oxidase (ChOx), a member of the glucose-methanol-choline (GMC) family, catalyzes the oxidation of the substrate via a hydride transfer mechanism and concomitant reduction of the FAD cofactor. Unlike other GMC enzymes, the conserved His447 is not the catalytic base that deprotonates the substrate in ChOx. Our QM/MM MD simulations indicate that the Glu361 residue acts as a catalytic base facilitating the hydride transfer from the substrate to the cofactor. We find that two rationally chosen point mutations (His447Gln and His447Asn) cause notable decreases in the catalytic activity. The binding free energy calculations show that the Glu361 and His447 residues are important in substrate binding. We also performed high-level double-hybrid density functional theory simulations using small model systems, which support the QM/MM MD results. Our work provides a basis for unraveling the substrate oxidation mechanism in GMC enzymes in which the conserved histidine does not act as a base.

Introduction

Cholesterol oxidase (ChOx) is a flavoenzyme that catalyzes the oxidation of cholesterol to 5-cholesten-3-one and the subsequent isomerization to form the final 4-cholesten-3-one product. The enzyme from Streptomyces is a member of the glucose-methanol-choline (GMC) family and contains a single flavin adenine dinucleotide (FAD) cofactor that is non-covalently bound to the protein. Apart from its physiological functions in bacterial metabolism, pathogenesis, and macrolide biosynthesis1, ChOx has also been employed as a useful biotechnological tool, for example, for the determination of serum cholesterol levels2,3.

A number of experimental studies investigated the structure and catalytic mechanism of ChOx over the past two decades4,5,6,7,8,9,10,11. The active site of ChOx is shown in Fig. 1a. With regard to the enzymatic mechanism, Sampson et al. identified three residues (Glu361, His447, and Asn485) that are significant for the enzymatic activity and Vrielink et al.4,5 revealed that Gly120 also plays an important role in the catalytic mechanism.

Figure 1

Key residues in the active site of ChOx and proposed reaction mechanism. ( a) The substrate is colored in magenta. The grey dashed lines represent H-bond interactions. Atomic color scheme: H, white; C, green; O, red; N, blue. ( b) Proposed mechanism for the substrate oxidation reaction catalyzed by ChOx.

Mutations of the Glu361 residue in the active site to glutamine and aspartate result in 31-fold9 and 14-fold8 reductions, respectively, in the k cat value for the oxidation of the substrate compared to the wild-type (WT) enzyme. This suggests that the carboxylic acid residue at this position may be important for the oxidative activity6. The His447 residue is located in the active site and the NE2 hydrogen forms a hydrogen bond with the hydroxyl group of the substrate. Mutations of His447 to glutamine and asparagine result in 140-fold and 4400-fold reductions in the k cat value respectively9, relative to the WT enzyme, suggesting that His447 plays a significant role in the catalytic activity. However, the double mutant, His447Gln/Glu361Gln, exhibits a 500-fold decrease in the k cat value relative to the WT enzyme10, which was 3-fold slower than that for the His447Gln single-mutant. These mutagenesis results indicate that factors affecting the oxidation rate of the two separate mutants, His447Gln and Glu361Gln, are not entirely additive and may be the result of structural perturbations rather than the absence of an active site base. We note that the double mutation does not completely prevent the hydride transfer implying that possibly the FAD is sufficiently electrophilic to oxidize the substrate without complete proton abstraction6. In addition, other structural studies have suggested that Asn485 forms an N-H***π electrostatic interaction with the flavin π-system to facilitate reduction of the FAD11. Mutation of Asn485 to leucine decreases the k cat value as well. It has been proposed that this residue modulates the electrostatic potential of the flavin to enhance substrate oxidation11. The Gly120 residue, which is positioned below the isoalloxazine ring system and functions as a hydrogen bond donor to the N5 atom of FAD, was also proposed to be important in orienting the orbitals of the N5 atom of the FAD cofactor, thus priming it for the hydride transfer chemistry4,5.

The recent high resolution X-ray structure of the oxidized ChOx in the presence of an isopropanol substrate (PDB ID: 4U2T)4 enables us to build a model of the ChOx active site composed of Glu361, the FAD cofactor, and the steroid substrate (see Fig. 1b). It was proposed that ChOx catalyzes the oxidation of the substrate via a hydride transfer from the steroid substrate C3 atom to the N5 atom of FAD and concomitant reduction of the FAD cofactor5,6. Despite extensive experimental work on the structure of ChOx, a detailed theoretical simulation of the catalytic mechanism and roles of the key residues in the active site was not carried out. In this work, we performed extensive combined quantum mechanics and molecular mechanics molecular dynamics (QM/MM MD) simulations including the catalytic mechanism and substrate binding, and high-level double-hybrid density functional theory (DHDFT) calculaions to explore the enzymatic catalysis. The results from the DHDFT calculations using small model systems support the QM/MM MD results. The combination of these two approaches gives an in-depth understanding of flavoenzyme catalysis by ChOx.

Results

QM/MM MD simulations on enzymatic mechanism

Catalytic mechanism of ChOx

10 ns molecular dynamics (MD) simulations for the enzyme-substrate complex model were performed. The root-mean-square deviation (RMSD) of the protein backbone is indicated in Supplementary Fig. S1. The last snapshot of MD simulations was chosen as the initial QM/MM model by removing the water molecules beyond 30 Å from the N9 atom of the FAD cofactor (see Fig. 2). We considered two different reaction coordinates (RC1 and RC2, see Supplementary Fig. S2). The reaction coordinate RC1 is defined as the difference between the C3-H3 and N5-H3 bond distances that are involved in the hydride transfer: RC1 = d C3-H3 - d N5-H3. The reaction coordinate RC2 is defined as RC2 = d C3-H3 - d N5-H3 - d O2-H1, that is it describes the movement of both the hydride and the proton. The relative energy barrier heights for RC1 and RC2 are 130.6 and 125.4 kJ mol −1, respectively. Another way of describing the combined movement of both the hydride and the proton is to fix the length of the O1-H1 distance at certain values and to scan along the reaction coordinate RC1. In particular, the O1-H1 distance was fixed at 1.01, 1.06, 1.11, 1.16, 1.21, 1.26, 1.31, 1.36, 1.41, 1.46, 1.51, 1.56, 1.61, 1.66, 1.71, and 1.76 Å. This is a cost-effective way to approximate a 2D scan of both the hydride and proton transfers. This approximated 2D scan is shown in the Supplementary Fig. S3a. It clearly shows that a first-order saddle point is obtained for the hydride transfer for each of the fixed O1-H1 bond distances. We find that the lowest reaction barrier height of 128.4 kJ mol −1 is obtained at r(O1-H1) ≈ 1.56 Å (Supplementary Fig. S3b). This result is consistent with our results for RC1 and RC2 above. Considering the computational cost of the potential of mean force (PMF) calculations, we chose RC1 as our final reaction coordinate since the reaction barrier height obtained along this coordinate is in good agreement with our best reaction barrier height obtained from the approximate 2D scan. By using QM/MM MD simulations and umbrella sampling, the free energy profile and corresponding configurations of reactant, transition state, and product involved in the enzymatic catalysis have been determined. In order to verify the efficiency of our sampling, the configurational overlap in the range of 5 ps to 20 ps is shown in Supplementary Fig. S4. This figure shows that the sampling covers the whole reaction coordinate with the current step size and biasing harmonic potential.

Figure 2

Overall view of the QM/MM model for ChOx and the active site. The active site includes the Lys225, Asn119, Gly120, Glu361, Asn485, and His447 residues as well as FAD and the substrate. The residues in the QM region are shown as ball and stick, the protein as ribbons, and the solvent as wireframes. The QM/MM boundary settings are specified in the subsection of "QM/MM MD simulation".

Figure 3 shows the free energy profile and the structures of the reactant ( React), transition structure ( TS), and product ( Prod) obtained from the umbrella sampling technique. It is noteworthy that the deprotonation of the substrate by Glu361 occurs synchronously with the hydride transfer (see TS in Fig. 3). The O1-H1 bond distance increases from 1.01 ± 0.03 ( React) to 1.90 ± 0.17 Å ( Prod), while the O2-H1 bond distance decreases from 2.10 ± 0.45 ( React) to 0.99 ± 0.01 Å ( Prod). The free energy barrier of the hydride transfer is 99.3 kJ mol −1, and the reaction energy is exergonic by 11.8 kJ mol −1. In the QM/MM transition structure, the length of the C3-H3 bond that is being broken is 1.37 ± 0.07 Å and the length of the N5-H3 bond that is being formed is 1.35 ± 0.09 Å. In this transition structure the O1-H1 and O2-H1 bond distances are 2.04 ± 0.30 and 1.00 ± 0.03 Å, respectively. This indicates that the proton has transferred before the hydride transfer has been completed. We note that we cannot directly compare the barrier obtained from the QM/MM MD simulations with previous experimental kinetic data9 since the kinetic measurements do not reflect the rate of the hydride transfer alone. However, the experimental k cat value obtained for the wild type (44 ± 2 s −1) indicates that our QM/MM MD free energy barrier for the hydride transfer likely represents an overestimation. Nevertheless, it is important to stress that the present work is concerned largely with relative reaction barrier heights (that will benefit from some cancellation of errors) and are likely to be more accurate than the absolute reaction barrier heights.

Figure 3

Structures of reactant ( React), transition state ( TS), product ( Prod), and free energy profile of the hydride transfer in ChOx. Hydrogen bonds are shown as grey dashed lines. The bonds being broken and formed are shown as black dashed lines. The average distances (in Å) for selected bonds are taken from umbrella sampling.

In order to clarify the order of the proton and hydride transfers in the catalytic process, the entire snapshots along RC1 have been collected. Figure 4 plots the snapshots for the C3-H3, N5-H3, O1-H1, and O2-H1 distances along the reaction coordinate RC1 (between the interval RC1 = −1.59 to 2.61 Å). As shown in Fig. 4a, the length of the C3-H3 distance increases smoothly from ~1.1 to 3.8 Å, while the length of the N5-H3 decreases from ~2.9 to 1.0 Å along RC1. The cross-point of these two lines corresponds to the TS shown in Fig. 3. At this point, the C3-H3 and N5-H3 distances are equal to one another (i.e., RC1 = 0.0 Å). Figure 4b shows that the proton transfer from the hydroxyl group of the substrate to the oxygen (O2) of Glu361 occurs at RC1 ≈ −0.2 Å. More specifically, the O1-H1 bond length remains at ~1.0 Å and the O2***H1 distance fluctuates at ~1.8 Å along RC1 from around RC1 = −1.7 to −0.2 Å. The atom H1 jumps to O2 at RC1 ≈ −0.2 Å, after which the O2-H1 bond length stays at ~1.0 Å and the O1***H1 length fluctuates at ~ 2.3 Å. These results indicate that the proton transfer occurs spontaneously with the hydride transfer.

Figure 4

Statistics for bond lengths from the snapshots along RC1. ( a) the C3-H3 and N5-H3 ( b) the O1-H1 and O2-H1 bond lengths.

The bond lengths of His447NE2H-O1 and LysNH ₃-O were also monitored during the hydride transfer (see Supplementary Fig. S5). The hydrogen bond distance between the NH ₃ group of Lys225 and the O=C group of Asn119 (labeled as LysNH ₃-O in Supplementary Fig. S5) fluctuates at ~1.8 Å during the hydride transfer, which implies that this hydrogen bond has a negligible effect on the barrier height for the hydride transfer. As for the hydrogen bond distance between His447 and the hydroxyl group of the substrate, labeled as NE2H***O1, it is ~1.8 Å before the proton transfer (RC1 = -0.2 Å) and increases to ~2.6 Å after the proton transfer. This indicates that the His447 plays a significant role during the oxidation reaction. The hydrogen bond formed between His447 and the hydroxyl group of the substrate orients the hydroxyl H atom towards the Glu361 base, which facilitates the proton transfer. This also results in an energetically favorable trans arrangement between the substrate hydroxyl H atom and the hydride.

The effects of mutating His447 on the active site environment

Previous experimental point mutations revealed that the His447 residue plays an important role in the oxidation reaction of the substrate in ChOx9. For example, mutations of the His447 residue to glutamine and asparagine result in 140-fold and 4400-fold reductions, respectively, in the k cat value for the oxidation reaction9. However, the atomic resolution crystal structure showed that the NE2 atom of His447 was protonated6. This led to the suggestion that His447 is not a catalytic base. To further elucidate the role of His447, QM/MM calculations have been performed on systems in which the His447 residue has been mutated to glutamine and asparagine (denoted by His447Gln and His447Asn, respectively). In these computational experiments, we carried out QM/MM energy scans along the reaction coordinate of RC1 using the QM/MM optimized structures to obtain the minimum energy reaction path (see Supplementary Fig. S6). However for reasons of computational cost, the subsequent QM/MM free energy perturbation and 20 ps QM/MM MD simulations with umbrella sampling were not performed in this section. Figure 5 gives the TSs of the WT, His447Gln, and His447Asn systems. The energy barrier height for the hydride transfer in the mutated systems is higher by 32.8 (His447Gln) and 66.3 kJ mol −1 (His447Asn) relative to that in the WT system (as shown in Supplementary Fig. S6). These QM/MM scan results are consistent with the mutagenesis experiments, where mutation of His447 to asparagine results in a larger reduction in the k cat value compared to the glutamine mutant.

Figure 5

Transition structures for the WT, His447Gln, and His447Asn systems obtained by means of QM/MM scans. Hydrogen bonds are shown as grey dashed lines. The bonds being broken and formed are shown as black dashed lines. The distances for selected bonds are given in Å.

The selected bond distances for the TSs of the WT, His447Gln, and His447Asn systems are also labeled in Fig. 5. In the QM/MM-optimized structures the hydrogen bond distances between O1 of the substrate and NE2H of His447, His447Gln, and His447Asn are 1.76 (WT), 1.87 (His447Gln), and 3.72 Å (His447Asn). As we discussed earlier, the hydrogen bond formed between His447 and the hydroxyl group of the substrate facilitates the proton transfer. When the His447 residue is mutated to asparagine, due to the shorter side chain in asparagine relative to histidine, a hydrogen bond to the substrate is not formed (i.e., the distance between O1 and NE2H is 3.72 Å). These results indicate that the weak hydrogen bond or lack of hydrogen bond interactions between the substrate and the key residue (i.e. His447) may impede the proton transfer that facilitates the hydride transfer in the WT enzyme, and further increases the barrier height for the hydride transfer. This could be due to the fact that the OH group of the substrate is not as well fixed in the active site and the catalytic environment has changed in the mutagenesis systems (i.e. His447Gln and His447Asn).

In order to consider the protein/water environment for the His447Gln mutant system, we carried out 20 ns MD and QM/MM MD simulations followed by a potential energy scan. Firstly, 20 ns MD simulations were performed for the His447Gln mutant and a stable mutant was obtained. The RMSD of the protein backbone is shown in Supplementary Fig. S8a. Then the QM region was chosen similarly to that of the wild-type and optimized at B3LYP/6-31 G(d) level of theory. In order to obtain a stable protein/water system, before running the potential energy surface scan, we ran 5 ps QM/MM MD for the His447Gln mutant system. Then we carried out the relative energy scans using the different reaction coordinates (RC1 and RC2, vide supra). The scan plots are displayed in Supplementary Fig. S9. The relative energy barrier heights obtained from RC1 and RC2 are 158.7 and 162.8 kJ mol −1 respectively, which are in good agreement with the one shown in Supplementary Fig. S6. We believe that these relative energy barrier heights demonstrate the important role that His447 plays in the enzyme. We note that the main purpose of this work is to show that the predicted theoretical trend of the barrier height when mutating His447 to Gln is consistent with the experimental reduction of the k cat value, rather than to reproduce the experimental k cat value, which would require extensive additional PMF calculations.

Similar to our approach with the wild-type, we confirm that the chosen reaction coordinates (RC1 and RC2) give results that are consistent with the results of the more rigorous approximated 2D scan for the His447Gln system. In particular, the O1-H1 distance was fixed at 0.98, 1.03, 1.08, 1.13, 1.18, 1.23, 1.28, 1.33, 1.38, 1.43, 1.48, 1.53, 1.58, 1.63, 1.68, 1.73, 1.78, 1.83, 1.88, 1.93, 1.98, 2.03, 2.08, 2.13, and 2.18 Å for each scan of the hydride transfer reaction coordinate. The 2D scan is given as Supplementary Information (Supplementary Fig. S11a). This figure shows that a first-order saddle point is obtained for the hydride transfer for each of the fixed O1-H1 bond distances. We find that the lowest reaction barrier height of 156.8 kJ mol −1 is obtained at r(O1-H1) ≈ 1.63 Å (Supplementary Fig. S11b). This result is consistent with our results for RC1 and RC2 above. Thus, we conclude that both the RC1 and RC2 reaction coordinates can efficiently describe the hydride transfer as well as the proton transfer. Considering the computational cost of the PMF calculations and the main purpose of this work, we chose RC1 as our final reaction coordinate for the potential energy scan.

The role of key residues in substrate binding

It has been proposed that the hydrogen-bonding networks involving Gly120, Glu361, His447, Asn485, and the substrate play important roles in substrate binding. It is interesting to evaluate the binding free energy and the interaction between the substrate and the key residues. The ligand-residue interaction for residues from 101 to 499 in ChOx and binding free energy decomposition in the active site are displayed in Supplementary Fig. S7. We note that Gly120, Glu361, His447, and Asn485 play important roles in stabilizing the substrate in the active site (see Supplementary Fig. S7a). It is noteworthy that other residues also play significant roles in the binding free energy, e.g., Met114, Met122, Val250, Pro364, Tyr446, Pro448, and Phe487. For example, they are associated with binding free energies ranging from −5.2 (Pro364) to −19.0 kJ mol −1 (Pro448). These residues are involved in the surrounding hydrogen bond network of the FAD cofactor and stabilize FAD in the active site. However, they are not directly associated with stabilizing the substrate in the active site.

In order to get additional insights into individual contributions from the Glu361, His447, and Asn485 residues to the substrate binding, the binding free energies were calculated for the WT and the Glu361Ala, His447Ala, and Asn485Ala mutant systems using the MM-GBSA method13,14. The Gibbs free binding energies for the WT and mutant systems as well as the differences between them are listed in Table 1. The Gibbs free binding energies for the mutants are lower than that for the WT by: 35.8 (Glu361Ala), 41.3 (His447Ala), and 4.8 kJ mol −1 (Asn485Ala). Mutations of these key residues result in destabilization of the substrate in the reaction center. Thus, the catalytic process becomes less efficient. This is another indication of the significant roles of these key residues in stabilizing the substrate in the active site. In order to obtain further insights into the roles of key residues in the enzymatic catalysis we also carried out double-hybrid density functional theory (DFT) calculations, which will be discussed in the following section.

Table 1: Gibbs binding free energies (∆G bind, in kJ mol−1) for the WT, Glu361Ala, His447Ala, and Asn485Ala systems as well as the differences between the WT and mutant systems (∆∆G bind) calculated with the MM-GBSA method.

Double-hybrid DFT calculations for probing the effects of hydrogen-bonding interactions on the barrier height of the hydride transfer

Double-hybrid density functional theory (DHDFT) calculations were performed in order to obtain insights into the roles of the residues in the catalytic mechanism of the hydride transfer from the C3 atom of the substrate to the N5 atom of the FAD cofactor. These computational experiments were carried out on a number of model systems with increasing size in order to elucidate the catalytic roles of the Glu361 and Gly120 residues on the barrier height for the hydride transfer. In particular, we consider the following models. Model A consists of the FAD cofactor (modeled by an isoalloxazine ring system) and the steroid substrate (modeled by an HOCHMe ₂ moiety). Model B consists of an additional Glu361, which forms a strong hydrogen bond with the hydroxyl group of the steroid substrate. Model C includes an additional Gly120 (modeled by a dimethylamine molecule), which forms a hydrogen bond to the N5 atom of FAD. The atoms involved in these three models are displayed in Supplementary Fig. S12.

Figure 6 depicts the Gibbs free energy profile (∆ G ₂₉₈) for the hydride transfer in models A (red), B (black), and C (blue). The calculated reaction barriers (∆ H ₂₉₈‡) and enthalpies (∆ H ₂₉₈) are listed in Supplementary Table S1. Figure 7 shows the optimized transition structures (TSs) involved in models A- C. The corresponding reactant complexes (RCs) and product complexes (PCs) are shown in Supplementary Fig. S13.

Figure 6

Gibbs free energy profiles for the hydride transfer in models A-C calculated at the RI-B2GP-PLYP-D3/Def2-QZVPP level of theory. Model A includes FAD and the substrate (red line); model B includes FAD, the substrate, and Glu361 (black line); and model C includes FAD, the substrate, Glu361, and Gly120 (blue line). The transition structures (TSs) involved in these reactions are shown in Fig. 7. The corresponding reactant complexes (RCs) and product complexes (PCs) are shown in Supplementary Fig. S13.

Figure 7

B3LYP-D3/6-31+G(2df,p) optimized TSs located on the potential energy profiles in models A, B, and C. Hydrogen bonds are shown as grey dashed lines. The bonds that are broken and formed in the TSs are shown in black dashed lines with purple circles. The C3-H3 and H3-N5 bond distances are given in Å. Atomic color scheme: H, white; C, grey; N, blue; O, red.

Model A, which consists of the FAD cofactor and the substrate, is the simplest binary system for modeling the hydride transfer. The hydride transfer in this model is expected to be both kinetically and thermodynamically unfavorable, however we start with this simple picture since it provides a reference point for the effects from the additional residues in models B and C. For the hydride transfer in model A, we obtain an activation Gibbs free energy of ∆ G ₂₉₈‡ = 136.5 kJ mol −1 and the reaction is endergonic by 37.1 kJ mol −1. According to the Hammond-Leffler postulate the high reaction barrier and energy suggest a late TS15. Indeed, in the TS, the hydride transfer between the C3 and N5 atoms is almost complete, the C3***H − and −H***N5 distances being 1.772 and 1.067 Å, respectively. For comparison, the −H***N5 bond length in the product complex is 1.035 Å (see Supplementary Fig. S13).

The enzymatic role of Glu361 is to act as a base to abstract the proton from the hydroxyl group of the substrate. Abstracting the hydroxyl proton has two effects on the reaction profile for the hydride transfer: (i) it activates the hydride transfer by pushing electron density to the C3-H3 bond and (ii) it stabilizes the product by converting the isopropanol carbocation to a neutral acetone. Model B includes the Glu361 residue. Indeed, the activation free energy in this model is lower by 22.6 kJ mol −1 relative to that in model A. In the TS of model B, the hydride transfer from C3 to N5 is concomitant with the proton transfer from the substrate to Glu361 (see model B in Fig. 7). In this TS, the hydride is located halfway between the C3 and N5 centers, in particular we obtain bond distances of 1.313 (C3***H ⁻) and 1.343 Å ( ⁻H***N5). These bond distances indicate that the TS for the hydride transfer in model B is earlier than that in model A, and therefore the lower barrier is consistent with the Hammond-Leffler postulate15. The effect of Glu361 on the barrier height of the hydride transfer is also demonstrated by considering the atomic polar tensor (APT) charges on the hydride and the substrate in the TS16,17. In model A, the overall charge on the substrate is +0.245 a.u. and the atomic charge on the hydride is −0.178 a.u. However, in model B, due to the proton transfer from the substrate to Glu361, the overall charge on the substrate becomes negative (namely it becomes -0.331 a.u.) and the hydride becomes more negatively charged as well with an atomic charge of -0.507 a.u. The partial negative charges on the hydride and the substrate in model B should facilitate an easier hydride transfer compared with model A where they have opposite charges. Abstracting the hydroxyl hydrogen in model B also stabilizes the product by converting the isopropanol carbocation to a neutral acetone. In this situation, the reaction energy is nearly thermo neutral, in particular it is reduced from +37.1 (model A) to +2.2 kJ mol −1 (model B).

It has been previously suggested that the hydrogen bond between the main chain NH of Gly120 and the flavin redox center (N5) may facilitate the catalytic hydride transfer by re-arrangement of the lone-pair electrons on N54. Model C includes a dimethylamine (HNMe ₂) molecule that mimics the hydrogen bond between Gly120 and FAD. The optimized TS is shown in Fig. 7 and the optimized O2***H1 and O1***H1 distances in the RC, TS, and PC are listed in Supplementary Table S2. The Gibbs free energy barrier height of model C is 102.1 kJ mol −1, i.e., it is 11.8 kJ mol −1 lower than that of model B (see Fig. 6), and the reaction energy of this model is exergonic by 14.7 kJ mol −1. These values are in excellent agreement with the barrier obtained from the QM/MM MD simulations, namely: ∆ G ₂₉₈‡ = 99.3 and ∆ G ₂₉₈ = -11.8 kJ mol −1. A natural bond orbital (NBO) analysis18 performed at the B3LYP-D3/6-31+G(2df,p) level of theory revealed that lone-pair → σ N-H interactions between the N5 center and N-H bond of Gly120 are associated with a stabilization energy of 20.6 kJ mol −1 in the reactant complex and 26.2 kJ mol −1 in the transition structure. The greater stabilization energy in the TS relative to the RC may partially account for the catalytic enhancement provided by the Gly120 residue. The above stabilization energies are approximated by the second-order perturbation energy ( E(2)) obtained from the NBO analysis19.

In summary, our high-level DHDFT modeling suggests that the Glu361 residue, which abstracts the proton from the hydroxyl group of the substrate, has a significant influence on the Gibbs free energy barrier for the enzymatic hydride transfer. Namely, the Glu361 residue lowers the barrier by 22.6 kJ mol −1 (comparison of models A and B). We also show that the Gly120 residue lowers the barrier further by 11.8 kJ mol −1 (comparison of models B and C).

Discussion

Extensive QM/MM MD simulations and DHDFT calculations have been carried out in order to explore the mechanism of the oxidation reaction catalyzed by ChOx as well as the roles that the surrounding residues play for the hydride transfer from the substrate to the co-factor.

Using QM/MM MD simulations and umbrella sampling we explored the reaction mechanism for the hydride transfer. In addition, the roles of the key residues in the substrate binding are examined. Point mutations of His447 to glutamine and asparagine increase the relative energy barrier height of the hydride transfer by 32.8 and 66.3 kJ mol −1, respectively. These results are in good quantitative agreement with the reduction in the experimental k cat values upon mutation. Our simulations suggest that mutation of His447 disrupts the hydrogen-bonding network around the substrate, which in turn results in the substrate being in a less favorable position for the proton transfer from the substrate to the base Glu361. In addition, the MM-GBSA calculations show that mutations of Glu361, His447, and Asn485 residues to alanine decrease the binding free energies of these systems by 35.8, 41.3, and 4.8 kJ mol −1, respectively, relative to the WT system. This indicates that Glu361 and His447 play important roles in substrate binding. On the basis of per-residue type, the free energy decomposition shows that the residues of Gly120, Glu361, His447, and Asn485 in the active site also have notable contributions to the substrate binding, in which the electrostatic and van der Waals interactions play predominant roles in binding the substrate in the active site.

The DHDFT results show that Glu361 and Gly120 play significant roles in the enzymatic catalysis. Glu361 serves as a base for abstracting the proton from the hydroxyl group of the substrate and facilitates the hydride transfer. The Gibbs free energy barrier for the hydride transfer is lowered by 22.6 kJ mol −1 upon inclusion of the Glu361 residue in our model. The hydrogen bond interaction between Gly120N-H and N5 of FAD also enhances the hydride transfer and lowers the barrier height of the hydride transfer further. We note that there is a good mechanistic agreement between the QM/MM MD simulations and our best DHDFT results.

In summary, our QM/MM MD simulations and high-level DHDFT using small model systems offer two complementary perspectives on the catalytic mechanism for the hydride transfer in ChOx. The combination of these two approaches gives an in-depth understanding of this flavoenzyme catalysis. More specifically, the QM/MM MD simulations show that the His447 and Glu361 residues play important roles in substrate binding. Our DHDFT models efficiently demonstrate the effects of the Glu361 and Gly120 residues on the reaction barrier height for the hydride transfer. The good agreement between the two approaches for the reaction barrier height of the hydride transfer increases our confidence in the results obtained from these models.

Methods

Preparation of the enzyme-substrate complex model

The initial coordinates used to build the model for the present study were based on the X-ray crystal structure of oxidized cholesterol oxidase (PDB ID: 4U2T)4. The steroid substrate dehydroepiandrosterone (DHA) was modeled in this structure in the proposed Michaelis state, based on the structure of the DHA/enzyme complex20 and knowledge of the interactions between amino acid side chains and the substrate from previous crystallographic and mutagenesis results. Specifically the substrate hydroxyl group was positioned so as to make hydrogen bond contact to NE2 of His447 and the C3 hydrogen atom was positioned above the flavin N5 atom. The protonation states of the ionizable residues were determined at pH = 5.2 with the H++ program22 and the neighboring hydrogen bond networks. The pH value used here is taken from experiment4. The FAD cofactor and the substrate as well as the protein were described by the AMBER GAFF force field23 and the AMBER14SB force field24, respectively. The electrostatic potential (ESP) was calculated at the HF/6-31 G(d) level of theory using Gaussian 09 program suite25. The computed ESP of FAD and the substrate were fit using restrained electrostatic potential (RESP) charge26 fitting in antechamber. The whole system was solvated into an 88 × 84 × 103 Å rectangular box of TIP3P water27 with a 12 Å buffer distance between the box edge and the nearest solute atoms. The system was neutralized by adding Cl ⁻ ions. The protons were added automatically, and the topology parameters as well as initial coordinates were generated by the tleap AmberTool. After the energy minimization, the system was heated gradually from 0 to 300 K for 100 ps, and another 100 ps MD simulation was carried out to relax the system density to about 1.0 g cm ⁻³. Finally, the standard 10 ns MD simulation under isothermal-isobaric (NPT) ensemble (NTP = 1) was employed with an integration time step of 2 fs via using the periodic boundary condition. The Berendsen barostat was used for pressure-control. The cutoff value was set to 10 Å for van der Waals and electrostatic interaction calculations. Long-range electrostatic interactions were dealt with by the Particle Mesh Ewald (PME) method28,29. The Langevin thermostat (NTT = 3) was used to maintain the temperature at 300 K. All the bonds involving hydrogen atoms were constrained by using the SHAKE scheme30. The RMSD was used to estimate the stability of the backbone of the enzyme, and the last 2 ns trajectories were utilized for the protein-ligand interaction decomposition of the binding free energy for both the WT and mutant systems (Glu361Ala, His447Ala, and Asn485Ala). Hou et al. reported that the MM-GBSA method with the generalized Born model performs well for binding free energy estimations31,32. In the current study, the MMPBSA.py module33 was then applied for the post-processing in which 1000 snapshots from an ensemble of conformations are used to calculate the free energy. This method has been also used to analyze the binding free energy for protein-ligand interactions for the WT and mutant systems. The per-residue decomposition was applied to estimate the contribution from each residue to the total binding free energy. The whole MD simulations were accomplished by applying the AMBER 14 software34.

QM/MM MD simulations

Development and applications of combined QM/MM methods for enzymes have been extensively reviewed in the past35,36,37. In order to obtain the barrier height of the hydride transfer, QM/MM MD simulations have been carried out for this enzyme, considering the influence from the surrounding residues to the reaction site. The residues: Lys225, Asn119, Gly120, Glu361, His447, Asn485, the FAD cofactor, and the substrate are considered as the QM region (103 atoms in total, see Fig. 2). The QM/MM boundary settings were determined according to the pseudo-bond rule38. These pseudo-bonds, C ps( sp 3)-C( sp 3), C ps( sp 3)-C( sp ², carbonyl) and C ps( sp 3)-N( sp 3) with accurate parameterization are preferred to be used when cutting the QM and MM regions. The specified sets are listed below in detail: (1) Lys227: only the -NH3⁺ and the adjacent -CH ₂ groups included (C ps( sp 3)-C( sp 3) type); (2) Asn119: only the -C = O group forming peptide bond with -NH group of Gly120 included (C ps( sp 3)-C( sp ², carbonyl) type); (3) Gly120: the whole Gly120 residue included and the -NH group of Gly121 included (C ps( sp 3)-N( sp 3) type); (4) Glu361: only the -COO ⁻ and the adjacent -CH ₂ groups included (C ps( sp 3)-C( sp 3) type); (5) His447: the imidazole ring with one -CH ₂ group included (C ps( sp 3)-N( sp 3) type); (6) Asn485: the -NH ₂ and two adjacent -C = O and -CH ₂ groups included (C ps( sp 3)-C( sp 3) type); (7) FAD: without the long chain (C ps( sp 3)-C( sp 3) type); (8) DHA: only one six-membered ring with the -OH group was included as well as two adjacent short chains (C ps( sp 3)-C( sp ², carbonyl) and C ps( sp 3)-C( sp 3) types). All the atoms in QM region are treated using the B3LYP/6-31 G(d) level of theory and atoms in the MM region are described by the AMBER99SB force field39. This approach has been successfully applied to other enzymatic reactions in the past decade40,41. The TIP3P model was used for the solvent water molecules27. The QM/MM boundary was addressed by the pseudo-bond approach37,38 with improved parameters42. The spherical boundary condition was applied, in which all atoms within a radius of 25 Å from the spherical center atom were allowed to freely move. The 12 and 18 Å cutoffs were used for van der Waals and electrostatic interactions among MM atoms, respectively. Following the interactive minimization, the minimum reaction energy path was mapped out by using the reaction coordinate driving method (RCD). Then, a 200 ps MD simulation was used to equilibrate the MM part with the QM subsystem constrained. The obtained snapshot was used for the subsequent QM/MM MD simulations combined with the umbrella sampling43. A total of 43 umbrella windows were employed to cover the reaction coordinate from -1.59 to 2.61 Å. For each window, 20 ps QM/MM MD simulations have been carried out with a time step of 1 fs by Beeman algorithm44. Berendsen thermostat method45 was used to control the system temperature at 300 K. The first 5 ps simulations were used for the equilibration and then another 15 ps simulations were used to obtain the PMF by the weighted histogram analysis method (WHAM)46. The appropriate step size and harmonic force constant for each window were chosen to guarantee the sufficient sampling. All the QM/MM MD simulations were carried out by the modified Q-Chem 4.047 and Tinker 4.2 programs48.

High-level double-hybrid DFT simulations

High-level DHDFT calculations using the B2GP-PLYP functional49 were performed in order to obtain accurate reaction energies and barrier heights of the hydride transfer from the C3 atom of the substrate to the N5 atom of the FAD cofactor. DHDFT methods include non-local electron correlation from second-order Møller-Plesset (MP2) perturbation theory in addition to the regular ingredients of hybrid DFT50. These methods overcome limitations of traditional DFT methods and display excellent performance for challenging chemical problems. Due to the slow basis set convergence of the MP2-type correlation component55, the DHDFT calculations were carried out in conjunction with the quadruple-zeta Def2-QZVPP basis set56.

The B3LYP-D3/6-31+G(2df,p) level of theory was applied to fully optimize the geometries in our models A, B, and C (see Fig. 7). Based on these optimized geometries, single-point energies at RI-B2GP-PLYP/Def2-QZVPP level were carried out using the ORCA program57 (RI stands for resolution of identity approximation applied in the MP2 step)58,59. Zero-point vibrational energy (ZPVE), enthalpic, and entropic corrections have been obtained from the B3LYP-D3/6-31+G(2df,p) harmonic frequency calculations for converting the electronic RI-B2GP-PLYP/Def2-QZVPP reaction energies and barrier heights into Gibbs free reaction energies and barrier heights at 298 K. Empirical dispersion corrections were included using the Becke-Johnson damping potential60 as recommended by Grimme et al. (denoted by the suffix D3)61,62,63. The equilibrium structures in the unconstrained optimizations were verified to have all real harmonic frequencies and the transition structures were confirmed to have only one imaginary frequency. The connectivities of the transition structures were confirmed by performing intrinsic reaction coordinate (IRC) calculations64,65. The enzyme-like environment was simulated by a homogeneous polarizable continuum model with the dielectric constant of 4.066,67,68 using the conductor-like polarizable continuum model (CPCM) at the HF/6-31+G(d) level of theory in conjunction with UAHF atomic radii, as recommended by Takano and Houk69.

All the above geometry optimization, harmonic frequency, and enzyme-like environment solvent correction calculations were performed using the Gaussian 09 program suite25.

Additional information

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. 5.

2. 24.

3. 25.
Frisch, M. et al. Gaussian 09, Revision D. 01. Gaussian: Wallingford, CT, USA, 2009.

4. 28.
Darden, T., York, D. & Pedersen, L. Particle mesh Ewald: An N⋅ log (N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089-10092 (1993).

5. 34.
Case, D. et al. AMBER 14, University of California, San Francisco, USA, 2014.

6. 39.

7. 48.
Ponder, J. TINKER: Software Tools for Molecular Design, 4.2, Washington University School of Medicine, Saint Louis, MO (2004).

8. 60.
Becke, A. D. & Johnson, E. R. A density-functional model of the dispersion interaction. J. Chem. Phys. 123(154101), 1-9 (2005).

9. 61.
Grimme, S., Antony, J., Ehrlich, S. & Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 132(154104), 1-19 (2010).

Download references

Acknowledgements

We gratefully acknowledge the generous allocation of computing time from the National Computational Infrastructure (NCI) National Facility, and system administration support provided by the Faculty of Science at the University of Western Australia (UWA) to the Linux cluster of the Karton group, the provision of an Australian Postgraduate Award (to L.-J.Y.), and an Australian Research Council (ARC) Discovery Early Career Researcher Award (to A.K., project number: DE140100311). We also thank the National Science Foundation of China (to Y.Z., project number: 21603057), and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase).

Author information

Affiliations

1. School of Molecular Sciences, The University of Western Australia, Perth, WA, 6009, Australia

Li-Juan Yu

, Emily Golden

, Alice Vrielink

& Amir Karton

2. The Key Laboratory of Natural Medicine and Immuno-Engineering, Henan University, Kaifeng, 475004, China

3. School of Pharmaceutical Sciences, Sun Yat-sen University, Guangzhou, 510006, China

4. Department of Chemistry, University of California, Davis, California, 95616, United States

Contributions

Conceptualization: L.-J.Y., A.V. and A.K. Methodology: L.-J.Y., N.C., Y.Z. and A.K. Formal Analysis: L.-J.Y., N.C. and Y.Z. Investigation: L.-J.Y., N.C., Y.Z., A.V. and A.K. Writing Original Draft: L.-J.Y. Writing Review and Editing: L.-J.Y., E.G., N.C., Y.Z., A.V. and A.K. Figure Preparation: L.-J.Y. and E.G. Supervision: A.V. and A.K. Project administration: A.K.

Competing Interests

The authors declare that they have no competing interests.

Corresponding authors

Correspondence to Yuan Zhao or Amir Karton.

Electronic supplementary material

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

# data[].attributes.attachment_link.title
Computational insights for the hydride transfer and distinctive roles of key residues in cholesterol oxidase

# data[].attributes.attachment_link.description
Altmetric: 7 More detail Article | Open Scientific Reports 7, Article number: 17265 (2017) doi:10.1038/s41598-017-17503-x Download Citation Received: Accepted: Published online: Cholesterol oxidase (ChOx), a member of the glucose-methanol-choline (GMC) family, catalyzes the oxidation of the substrate via a hydride transfer mechanism and concomitant reduction of the FAD cofactor.

# data[].attributes.body
Geometry Optimization with Machine Trained Topological Atoms

François Zielinski, Peter I. Maxwell, Timothy L. Fletcher, Stuart J. Davie, Nicodemo Di Pasquale, Salvatore Cardamone, Matthew J. L. Mills & Paul L. A. Popelier
06 September 2017
doi:10.1038/s41598-017-12600-3

The geometry optimization of a water molecule with a novel type of energy function called FFLUX is presented to bypass the traditional bonded potentials. FFLUX’s architecture is suitable to predict their IQA atomic energies for a previously unseen molecular geometry. Training of FFLUX is achieved via an in-house script called GAIA, which controls the construction of a kriging model. The GAIA protocol has five key steps: (1) sampling, (2) ab initio calculations, (3) atomic property calculations, (4) kriging model building and (5) validation. The steps are carried out sequentially, with the output of the previous step forming the input for the next step.

# data[].attributes.attachment_link.readable_text
Altmetric: 1

Citations: 1

More detail

Article | Open

Scientific Reports 7, Article number: 12817 (2017)

doi:10.1038/s41598-017-12600-3

Download Citation

Received:

Accepted:

Published online:

Abstract

The geometry optimization of a water molecule with a novel type of energy function called FFLUX is presented, which bypasses the traditional bonded potentials. Instead, topologically-partitioned atomic energies are trained by the machine learning method kriging to predict their IQA atomic energies for a previously unseen molecular geometry. Proof-of-concept that FFLUX's architecture is suitable for geometry optimization is rigorously demonstrated. It is found that accurate kriging models can optimize 2000 distorted geometries to within 0.28 kJ mol −1 of the corresponding ab initio energy, and 50% of those to within 0.05 kJ mol −1. Kriging models are robust enough to optimize the molecular geometry to sub-noise accuracy, when two thirds of the geometric inputs are outside the training range of that model. Finally, the individual components of the potential energy are analyzed, and chemical intuition is reflected in the independent behavior of the three energy terms (intra-atomic), (electrostatic) and (exchange), in contrast to standard force fields.

Introduction

Traditional force fields express energy as a function of the internal coordinates of a chemical system. These potential energy functions are only loosely connected to an underlying quantum mechanical reality, if at all. Typically, the various force field energy contributions each fall into one of two broad categories: bonded (covalent) and non-bonded (non-covalent). Although this may appear a natural and innocent partitioning, the sharp distinction does not properly reflect the complexity of the atomic interactions found in condensed matter. Hydrogen bonding is probably the oldest type of interaction to challenge the artificial distinction between bonded and non-bonded interaction. Indeed, the modern approach of Interacting Quantum Atoms (IQA)1, which works with finite-volume topological atoms2,3,4, offers a view of covalency as a sliding scale5,6. Despite the built-in, artificial nature of their composite functions, popular force fields preserve the binary approach, and utilize a variety of bonded energy terms (such as bond-stretching, angle-bending, torsional rotation and their cross terms), and non-bonded energy terms (such as van der Waals interactions and point charge electrostatics).

A second major feature of these force fields is that the energy expressions are written as penalty functions. For example, if a given bond takes on its equilibrium bond length then the corresponding bond stretch energy is zero. Any deviations from equilibrium (either by bond compression or elongation) result in a positive energy penalty. The force field thus needs a reference geometry (i.e. the equilibrium geometry). Furthermore, the typical Lennard-Jones potential appearing in the modelling of van der Waals interactions introduces its own reference minimum-energy distances. The electrostatic interaction, which is typically written as a Coulomb interaction between point charges, introduces another reference, namely, that of charges being infinitely far apart.

Here, we use a very different approach, called FFLUX. This method, which was formerly called QCTFF7, is one in which atoms endowed with quantum mechanical knowledge8 come together to form a molecule. The topological energy partitioning method, IQA, offers a route to accomplish this goal when it is combined with a machine learning method. The latter (in this work) is kriging9,10,11 (or Gaussian process regression), which unlike neural networks or genetic algorithms originated in geostatistics. Kriging is a method of interpolation, giving the best linear unbiased prediction of the intermediate values. In 2009 this method was first used12 in combination with topological atoms in work that successfully captured the fluctuation of multipole moments of atoms in water clusters (up to the hexamer) in response to geometrical changes in the clusters. This advance constituted the first application of kriging in the context of intermolecular potentials, soon followed by the careful construction13 of interatomic potentials for solid state simulations.

We further developed the aforementioned proof-of-concept, and have demonstrated the applicability of kriged topological atoms in a growing variety of cases, including: water clusters14,15, methanol16, N-methylacetamide17, cholesterol18, a microhydrated sodium ion19, all proteinogenic amino acids20,21 including aromatic amino acids22, alanine helices23, hydrogen-bonded24 and weakly bound complexes25 (both from the S22 data set), and carbohydrates26. This collective work displayed the performance of the kriging models in terms of the accuracy of their energy predictions. For that purpose, we typically showed the cumulative error distribution (the so-called "S-curve") of the energy of each of the test geometries of the kriged system. As the intersection between the set of training geometries and the set of test geometries is the empty set, this type of validation is external. Over the years, the validation of the kriging models has been very systematic, complete and candid. In this same tradition we now systematically investigate a truly novel type of geometry optimization where a kriging model informs an atom on "how to behave" in the presence of other atoms. Ever since the availability of analytical forces27 for kriging models, it has been possible to make nuclei move towards an energy minimum. Herein the first application of this technology is described.

In the current paper we show how the geometry of a water molecule can be optimized without ever using bonded force field potentials. Instead, we let atoms adjust themselves as guided by quantum mechanical energy contributions as defined within IQA, including kinetic, exchange, and Coulomb energies. We focus on the essential, technical points of the method, and on its application. We then report various observations collected from a variety of kriging water models: (i) a statistical assessment of the kriging models' predictions, (ii) an in-depth look at the optimization performance of the various models and optimizers' parametrizations, (iii) a robustness test over a large set of starting points, and (iv) chemical insight obtained from the IQA framework in combination with FLUXX.

Computational Methods

There are several components of the FFLUX approach that cannot be discussed in great detail here because of space limitations. Below we provide key references to previous work where these components have been carefully and extensively explained. Training of FFLUX is achieved via a number of in-house and external computer programs, called by an in-house script called GAIA, which controls the construction of a kriging model. A detailed flowchart of GAIA is given in the appendix of reference22. GAIA controls thousands of input and output files and allows a user to essentially parameterize FFLUX for any system of interest. The GAIA protocol has five key steps: (1) sampling, (2) ab initio calculations, (3) atomic property calculations28, (4) kriging model building and (5) validation. Each step is carried out sequentially, with the output of the previous step forming the input for the next step. Before discussing these five steps in turn, we summarize in Fig. 1 the programs used in the construction of kriging models. The first four steps of the GAIA protocol correspond with the first four steps in the diagram of Fig. 1, while the validation step in GAIA usually takes the form of an S-curve (which establishes a kriging model's energy prediction quality) but can also take the form of an geometry optimization, and later in FFLUX's development, of an evaluation of a (thermo)dynamic property emerging from a condensed matter molecular dynamics simulation. In the final step of Fig. 1 we mention the molecular dynamics program DL_POLY, which can perform these validations, because it uses FFLUX's kriging models (in so-called production mode).

Figure 1

Flowchart of FFLUX's training (first four steps) and execution (DL_POLY), detailing the programs involved and summaries of their tasks.

Sampling by normal modes

A single water molecule was geometry optimized at HF/6-31 + G(d,p) using the program GAUSSIAN0929 (all default settings, with " nosymm" and " 6D"). A single-point frequency calculation was then performed to compute the second derivatives of the potential energy (Hessian) at the optimized geometry so that the normal modes of vibration could be determined. As described in Cardamone et al.26, the optimized geometry can be used as a "seed", about which one can stochastically sample the molecular normal modes. One is then able to collect an ensemble of such samples, and use them as the input for kriging, as we describe in a later section. This sampling methodology has been implemented in the in-house code EROS. To prevent any unphysical geometries from arising during the sampling process, constraints are placed on any samples included in the ensemble. We require that all bond lengths and valence angles are "distorted" by no more than ±20% of their values in the original "seed" geometry. So, for instance, if a given bond length is 1.0 Å in the "seed" geometry, we constrain the value of the bond length to lie in the range 0.8 Å-1.2 Å for all samples in the ensemble. While the magnitude of distortion is a free parameter, we have found from experience that a value of ±20% allows for an extensive sampling of conformational space, without having to resort to the multi-reference wavefunction techniques required for heavily distorted systems. This procedure was used to generate 2000 water geometries distorted from the geometry-optimized seed.

Ab initio calculations

The wavefunction of each of the 2000 geometries was calculated by GAUSSIAN09 at HF/6-31 + G(d,p). Here we are interested in demonstrating that FFLUX geometry optimization works. Hence we proceed with Hartree-Fock wavefunctions, whose IQA partitioning is simple and well-defined; more advanced wavefunctions will be introduced in subsequent work. Note that the QM minimum wavefunction is never included in the machine learning training sets described below, in order to more strictly test the capability of kriging.

Atomic property calculations

The IQA method is part of an overall approach, coined30 Quantum Chemical Topology (QCT)31 in 2003, which is based on the central idea of (gradient) vector field partitioning. This crucial idea lies at the heart of the Quantum Theory of Atom in Molecules (QTAIM)2,32, which was the first component of QCT. Recently, QCT has been didactically explained from various angles3,33,34. QTAIM (and hence QCT) defines a topological atom, which has a well-defined electronic kinetic energy35. This feature is important in the design of a force field with a deep connection to quantum mechanics. IQA partitions a molecule's energy, , into a sum of atomic energies, ∑A E IQA A , which in turn are composed of intra-atomic and inter-atomic energy components,

where A and B represent atoms, and the subscript denotes the type of energy contribution. This equation contains the four types of IQA energy contribution that are relevant to the current study: the overall atomic energy , the intra-atomic (or self) energy , the exchange energy V x A B, and finally the (classical) Coulomb energy . We now briefly explain these primary energy contributions.

The intra-atomic energy consists of kinetic, T, and potential energy, V, contributions:

where T A represents the kinetic energy of atom A, is the (attractive) potential energy between the electrons and nucleus of atom A, and is the (repulsive) potential energy between the electrons within atom A. The latter quantity can be generalized for any atom pair, , and further broken down as follows:

where 'Coul' refers to the Coulombic interaction between the electrons and 'x' represents the exchange energy. A third term, representing the correlation energy, is missing at Hartree-Fock level. Now that the Coulombic energy has been separated from , the classical electrostatic energy can be calculated by including the interaction involving the nucleus of A and of B,

Note that the order of the superscripts and subscripts is important because , for example, refers to the electrons of A interacting with the nucleus of B, and not the other way around.

Now that the explanation of the primary energy contributions is complete, one more remark needs to be made. A recent FFLUX publication36 introduced the use of interatomic energies designated by AA' instead of AB, a notation employed herein. Here A' represents every other atom in the molecular system except A. Thus, the notation AA' denotes the interatomic energy between an atom A and its surrounding environment A', such that

where the approximate equal sign is due to the AA' energy being calculated analytically (which is more accurate), but the AB energies by quadrature. Finally, the commercial package AIMAll37 is used to calculate these energy contributions from the wavefunctions using default parameters, except for the use of AIMAll's own implementation for computation of the two electron parts of (rather than the so-called TWOe integration option).

Kriging model building

Any machine learning method essentially links a set of inputs (called "features") with a set of outputs. Our use of kriging links a set of nuclear coordinates (the features or inputs) with a single output at a time (i.e. one of four possible types of atomic energy). Because the output depends only on the internal geometry of the molecule, there are 3 N-6 features, for a system with N atoms. In the case of water there will be three features: two O-H bond lengths and the H-O-H angle. The general definition of features introduces a broader context in which a local axis system, called an atomic local frame (ALF), must be installed on the atom being trained for. Strictly speaking, the axis system is only necessary for outputs that are directional quantities, such as atomic multipole moments, which do not appear in this article. However, the idea of installing an origin at the nuclear position of each atom, one at a time, must be explained here because this installation determines the way the features are constructed.

The x-axis of the ALF points from the origin atom to its heaviest bonded neighbour (following the Cahn-Ingold-Prelog convention). The xy-plane sweeps out from the x-axis toward the second heaviest atom bonded to the origin atom. The origin atom and first and second bonded atoms then determine the xy-plane. Subsequently, the y-axis is constructed to be orthogonal to the x-axis and the z-axis orthogonal to both, forming a right-handed axis system. The first three features consist of: (i) the distance between the origin atom (A ₁, 1 st atom) and the " x-axis-atom" (A ₂, 2 nd atom), (ii) the distance between the origin atom (A ₁) and the " xy-plane-atom" (A ₃, 3 rd atom), and (iii) the angle A ₂-A ₁-A ₃. For oxygen in water, the features are d(OH ₁), d(OH ₂) and α(H ₁OH ₂); for H ₁, they are d(OH ₁), d(H ₁H ₂) and α(OH ₁H ₂); and for H ₂, they are d(OH ₂), d(H ₁H ₂) and α(OH ₂H ₁). For large molecules, the molecular geometry is then converted from Cartesian coordinates in a global frame, to spherical polar coordinates of each remaining atom in the ALF (i.e. those atoms which are not part of the installation of the ALF). Note that each atom in the system acts as an origin for its own ALF, allowing the description of the remaining atoms by a unique (but complete) set of spherical polar coordinates.

Each atom in the system now sees its environment as a set of features (model inputs) and has a set of IQA properties (outputs, one per kriging model) that together make up a single training example for that atom. Since each geometry is a unique training example, and 2000 geometries were sampled, each atom in the system has a list of 2000 training examples (termed a 'sample set'). GAIA cleans the data through a "scrubbing" process by finding any examples with an AIMAll integration error (measured by L(Ω)) larger than a specified threshold, and removing these examples from the sample set. Any geometry that is incomplete due to removed atoms is then discarded from all atoms' sample set. Two atomic integration L(Ω)38,39 threshold values were used in this work. First, a value of 0.001 au was used for the training sets called "100", "300" and "500" in the following (which contain 100, 300 and 500 training examples, respectively). Second, a twenty times tighter threshold (of 0.00005 a.u. or 0.13 kJ mol −1) was applied for a second training set with 500 examples, termed "T500". Additionally, a third training set was created, termed "TE500", which replaces , , with a single value, namely that of . We have shown before that kriging can successfully36 construct a relationship between the various energy contributions and the (geometrical) features.

The kriging method outlined here is based on the treatment of Jones et al.40,41 and has been explained in much greater detail in our previous work17 and also in references14,42. Kriging maps the response of an output ŷ (an IQA energy) to any given input x (set of geometric coordinates),

where μ ˆ is the estimated global mean of the process, the background value for this output, and n is the number of training geometries. The quantity a i is the i th element of the vector a= R −1(y−1 µ ˆ) where R is a matrix of error correlations between training points, and 1 is a column vector of ones. The error from the global term is determined17 by the distance between the new input point ( x) and a known input point ( x i). The sum of these errors gives the appropriate deviation from the background term and results in the new output, y ˆ(x). An IQA energy has a mean 'background' value when considered across many geometries and kriging can map the deviations from the mean in response to geometric changes. The fact that kriging uses the distance between the new input and known inputs is chemically sensible as we can assume that if two geometries are very similar, the IQA energies on the atoms in each geometry are similar as well.

The symmetric correlation matrix R consists of the following kernel,

where d is the number of features, that is, the dimensionality of the input space, which is 3 in the current case study. In general, this value is equal to the number of internal coordinates, i.e. 3N-6. The correlation between two points in the training data is a function of the distance between the points, along with the kriging hyperparameters θ and p. These two sets of parameters may both be optimized in order for this correlation to best describe the effect that a move between these two inputs has on the selected output. Note that each dimension (feature) of the kriging problem has its own θ h and p h value. It has been suggested that p h can be fixed at 2 for most cases but instead it is optimized alongside θ h, because optimizing p h tends to help the kriging process with small molecules. This process is carried out with an implementation of the particle swarm optimization algorithm, with a log-likelihood objective function in our in-house code FEREBUS42. When the kriging training process is complete, a model is created that can be used to predict the IQA properties belonging to an atom when given a previously-unseen geometry. Any remaining data (that is, data that has not already been used to train models) left in the training sets can potentially be used as test examples as the kriging models have no knowledge of these examples.

Validation

The long-term strategy of FFLUX is one of bottom-up validation. This means that a kriging model is first assessed by the accuracy of its energy predictions. This is done via an S-curve, which will be explained in the Results Section. The next level of validation is based on prediction of a minimum energy geometry, which is achieved for the first time with FFLUX in this article. In order to achieve this aim, we need the analytical force that applies on each nucleus.. The next level of validation, which will be achieved in future work, is that of structure and dynamics obtained from a molecular dynamics (MD) simulation. This highest level of validation will appear in future work on kriged topological atoms, which covers the case of polarizable atoms if kriging trains for atomic dipole moments (and higher rank moments). In the case of non-polarizable topological atoms, their high-rank multipolar potentials have previously been tested against experiment for radial distribution functions and thermodynamic properties, for liquid water43,44, liquid imidazole45 and aqueous imidazole solutions46.

Atomic forces

For each Kriged quantity, i.e. IQA atomic energy, first derivatives are computed through adapted routines from earlier work dedicated to Kriged multipolar electrostatic interaction27. In the present case, in the absence of multipole moments, the reported mathematical framework essentially simplifies itself to:

where we differentiated with respect to the i th Cartesian coordinate α ( i = 1, 2, 3 referring to x, y, or z, respectively) expressed in the global frame. Each E and V quantity is represented by a kriging model or sum thereof. Note that a kriging model is expressed with respect to internal coordinates,

where E refers to any of the four types of energy, for any atom or pair of atoms, x = { x h; h = 1, 2, ..., d} is a given set of features or internal coordinates for which the energy needs to be predicted, and d and n again refer to the number of features and training examples, respectively. The derivative of this term with respect to an internal coordinate is given by:

where the cusp in the derivative of the absolute difference at = x k is dealt with by defining

The chain rule of differentiation serves as a bridge between Cartesian and internal coordinates or,

The kriging derivatives of the various energies can then be directly summed into atomic forces, once converted from the ALF to the global Cartesian frame by applying the chain rule.

Geometry Optimization

DL_POLY v4.0547 was chosen to host adapted code from the group's kriging prediction engine into a prototype dedicated to both proof-of-concept and design explorations. Several of the host software's capabilities, parallelism in particular, have been deactivated to facilitate the design and implementation of our method. The current local code demonstrates the viability of our method for gradient-based optimization techniques on water.

The current implementation chose to keep the new modules as self-contained as possible in order to minimize intervention into DL_POLY's core. By doing so, no changes had to be made to the Verlet integration or the optimization routines, which enabled seamless operation of the MD software. The 0 K (zero Kelvin) optimizer is equivalent to a molecular dynamics run set at minimal temperature, with the particles' velocities reset to zero between each step. In practice, an atom strictly moves along the forces to which it is subjected. Similar to a MD run, such an optimization then relies on the duration and number of timesteps as parameters. Long timesteps mean the optimum would be reached faster, at the risk of overshooting or oscillating around a narrow and deep minimum, while short timesteps would converge more slowly and risk being stuck in shallow and spurious local minima appearing in an undulating PES.

The conjugate gradient (CG) method proceeds by following the direction of a guess vector until the system's energy rises, at which point a new vector is computed, as a conjugate of the current last-guess. As a first parameter, the length of CG steps is based on the timestep length in DL_POLY's implementation. The latter also provides three different convergence criteria to stop the optimization process when satisfied: energy, forces, or displacement ("distance"). So far, only the last one is compatible with our kriging engine. Further details on the optimizers can be found in DL_POLY's user manual.

Results

Kriging model quality (S-curves)

The optimized geometry of water (hereafter referred to as the "QM minimum") was distorted using the program EROS with a maximal ±20% bond stretch and angle control parameter. The resultant O-H bond range was 0.754 Å ≤ x ≤ 1.132 Å, and the H-O-H angle range 85.70° ≤ x ≤ 128.54°. The sample set had a molecular energy range of 201.6 kJ mol −1. The program FEREBUS was used to obtain the molecular models. The quality of the five molecular models (100, 300, 500, T500 and TE500) is illustrated in the S-curves in Fig. 2, supplemented by the statistics given in Table 1. S-curves plot the prediction error ( x-axis) of each test prediction as a function of the number of test points (in our case 500, y-axis) so that each increment of 100%/500 = 0.2% on the y-axis represents a test point. Plotting test predictions on an S-curve allows a thorough inspection of a kriging model's quality. Within an S-curve, hallmarks of a good model are: (i) a steep gradient over a wide range centred at the curve's inflection point, (ii) being positioned as much as possible to the left, and (iii) a short 'tail' at the top (near 100%). The tail refers to the general portion of the curve where the highest errors are seen on the approach to the final point at 100%.

Figure 2

S-curves for the 100, 300, 500, T500 and TE500 water models described using the three energies given in Eqn. 1 (, and ). The label "T" stands for the tighter scrubbing threshold of 0.00005 Hartrees, while "TE" stands for this tight model using single total atomic energies, .

Three observations follow from the S-curves: (i) increasing the training set size from 100 to 300, and again to 500, incrementally moves the curve to the left, resulting in lower Mean Absolute Errors (MAEs), (ii) using a tighter scrubbing threshold T (0.00005 Hartrees instead of 0.001 Hartrees) showed little effect on the position of the S-curve, and (iii) kriging the single atomic energy instead of each of , and dramatically increased the gradient of the S-curve and shortened the tail. The statistics in Table 1 show us that model TE500 has a smaller range of errors, but a slightly higher MAE (0.10 kJ mol −1) compared to either 500 (0.06 kJ mol −1) or T500 (0.07 kJ mol −1). However, with the exception of the 100 model, all models performed very well, having very low molecular energy errors throughout.

The first section in the Supplementary Information shows an extensive cross-validation analysis showing the adequacy of the 300 water model to capture the behavior of the system.

Optimization Runs

Having ensured that the generated models are of good quality, the investigation now moves onto their application within the geometry optimization study. From the total sample set, three test samples were chosen as starting points (SP) for the initialization of DL_POLY's geometry optimization run. The three starting points are referred to as SP1, SP2 and SP3, and their relative molecular energies are, respectively, +15.05 kJ mol −1, +47.97 kJ mol −1 and +126.18 kJ mol −1 above that of the QM minimum. Selecting three starting points allows us to investigate each individually, but also to compare and contrast the resulting energies and geometries from each. The three SPs are selected to represent an incrementally more challenging task (from SP1 over SP2 to SP3). However, the geometries of each starting point (illustrated in Fig. 3) also feature three significantly different H-O-H angles of 115.62°, 106.3° and 93.36°, and also three quite different O-H bond length combinations. The increasing molecular energy of each starting geometry, along with three very different bond angles and bond length combinations (which notably feature no symmetry), ensure that the starting points begin their optimization trajectory from significantly different regions of conformational space. The described approach ensures a thorough assessment of each model.

Figure 3

SP1 (left, +15.05 kJ mol −1), SP2 (middle, +47.97 kJ mol −1) and SP3 (right, +126.18 kJ mol −1) water geometries. Bond distances are in Å, and bond angles in degrees.

For our preliminary exploration we select two different parameter sets for each of the two optimization methods. The 0 K optimization algorithm was run for 5000 steps with timesteps of 1 fs (parameter setting 1 or "Set 1" in short) and 0.5 fs ("Set 2"), while CG uses 1 fs timesteps with a distance convergence criterion of 10 −5 Å ("Set 3") and 10 −6 Å ("Set 4"). In total, four parameter sets are employed within this investigation

Finally, the kriged PES can be analyzed through a comparison of the geometric features of the optimized water molecules. Should the models have a similar kriged PES, the energies and geometries should show a similar optimization evolution, when starting from the same point. As a first step, here we will limit ourselves to quantitative comparisons on the optimized molecules only. Inconsistent results indicate that the respective PESs of the models are not so similar in the region of the optimal solution. From different SPs, any variation within the resulting energies and geometries within the same model will indicate an undulating PES(in which a trajectory can become trapped in a local minimum).

Table 2 summarizes the molecular optimization energy results for each SP, for all five molecular models. The QM energy (−199,620.00 kJ mol −1) is used as the reference energy for all Δ E values, calculated as: ΔE = [Final Molecular Energy −  QM Energy]. For each parameter setting (i.e. "Set"), the final optimized geometry energy is given along with the corresponding Δ E. We now discuss, in turn, four observations.

Table 2: SP1, SP2 and SP3 water E IQA Mol optimization results for each model (100, 300, 500, T500 and TE500).

The first observation is that all the optimizations have indeed run successfully with the exception of one (SP1, parameter set 4, for model 500 - details discussed shortly). However, before any further analysis, proof-of-concept has been shown: QCT atoms dressed up with IQA atomic energies and converted into kriging models, indeed are sufficient to obtain atomic forces suitable for molecular geometry optimization.

The second observation is that for all parameter sets (with only the above exception), the Δ E of the final geometry is ≤ ±1.6 kJ mol −1. In fact, most cases are ≤ ±0.2 kJ mol −1. Remarkably, the lowest Δ E reported is <0.01 kJ mol −1. However, final Δ E's smaller than the model's MAE (typically within ±0.1 kJ mol −1), are within the accuracy threshold of our approach. Also, the atomic integration implemented in AIMAll introduces energy noise that typically does not enable us to recover the QM energy closer than within ±0.1 kJ mol −1 (or 0.00005 a.u.). With the above in mind, it is remarkable to observe the kriging predictions performing so well. The same energetic minimum is being reached consistently for most SPs, within the same molecular model. Hence, no spurious local minima significantly corrupt the kriging PESs, which still appear unimodal. This observation can be made across models too but now returning more variation in the values of the energy minimum reached. Moving from the 100 model, over the 300 model to the 500 model, shows that denser sampling converges closer to the ab initio minimum. Note that there is pressure to keep the number of training examples to a minimum because of the computational cost of generating atomic energies. Hence, going beyond 500 cannot be justified, especially given the already excellent results obtained with 500. In summary, our aim is to build a kriging model with the smallest possible training set size, compatible with the accuracy we need for the problem considered.

We also see the abovementioned consistent convergence when operating both above and below the accuracy threshold. For example, in the TE500 model, all six 0 K optimizations (two across each of SP1, SP2 and SP3) all converge to a geometry with the same energy (Δ E = +0.14 kJ mol −1), above the accuracy limit (±0.1 kJ mol −1). For the 500 model we see this consistency (reaching Δ E = +0.04 kJ mol −1, below the accuracy threshold of ±0.06 kJ mol −1) for five of the six 0 K parameter settings. The respective geometries reached for the TE500 and 500 models, are not only energetically the same, but also geometrically, as shown in Table 3, which reports the final geometries of parameter Set 1.

The third observation is that the 0 K optimizations consistently perform better than the CG optimizations, where "better" equates to a lower Δ E. Throughout all 15 datasets (5 models × 3 SPs), there is only one example where both CG optimizations perform better than the two 0 K models (100 Model, SP3). The superior accuracy of 0 K was expected, given this method's design to follow the trajectory of a very low temperature simulation continuously until a chosen number of steps are reached. However, CG places more emphasis on reaching the minimum of the model more quickly, through a reduced number of steps and stopping when a convergence criterion is met. Indeed, although the use of CG has reached a slightly poorer molecular energy minimum geometry than 0 K, is the results are still good. Within the optimization community, the reduced accuracy of CG is generally accepted in favour of a fast calculation. Hence, CG's firm grounding as a common optimization algorithm, and 0 K's relative obscurity. For our investigation, 0 K is an undoubtedly useful algorithm for diagnostic purposes and proof-of-concept.

The fourth observation is that "Set 4" fails for SP1 for the "500" training as a result of the trajectory of the optimization never meeting the distance convergence criteria. Such behavior is an indication that the PES modelled in the IQA model is not smooth enough to reach a solution. At some point, the energy gradient of the model may cause a step in an incorrect direction of conformational space. Should such an event lead to the geometry 'escaping' the training range far enough to result in the kriging training correlation vanishing, then the geometry is considered to be in the flat "no man's land" that exists outside of the PES. For the SP1-Set 4 example, the molecular energy fluctuates between a good prediction and a poor prediction as the optimizer attempts to improve it. Eventually the model predicted a point far outside the training range, from which the trajectory failed to recover. A sensible hypothesis is that the PES is not accurate enough within that specific region of the conformational space. The lack of accuracy in the PES is likely either due to nearby points with high integration errors; to the fact that the training set used is not able to correctly describe the PES; or a mix of these two. If the problem is identified in the training set, the poor description of the PES could be due to an insufficient number of training points or a poor sampling of the conformational space. It is very difficult to tell a priori which is the main source of error and the more effective ways used to reduce this problem include (i) the use of more training points, (ii) considering more accurate training points, (iii) a better distribution of training points in conformational space, (iv) a combination of the three solutions presented before. Alternatively, completing the optimization with a less strict distance convergence criterion (as in "Set 3") is a solution that would not involve the modification of the model.

Having analyzed the energy of each final timestep, we now analyse the energy evolution trajectory. Because we observe consistently low energy Δ E values for the T500 model, these energy trajectories should provide a good example of the behaviour that one can expect from an accurate kriging water model. Figure 4 shows the energy trajectory of the optimization using Set 1 (0 K and 1 fs) and Set 3 (CG and 10 −5 Å convergence) for each starting point. The left panels show a difference between the smooth trajectory of the 0 K algorithm (top) and the jagged trajectory expected from CG (bottom). This behaviour is amplified in the right-hand magnified plots. The magnified plots show the convergence for each SP through monitoring the energy differences between successive timesteps. The plots on the right of Fig. 4 define Δ E as [ current system energy − previous system energy]. Note that the convergence plot for CG always ends with a peak, caused by the CG algorithm being forced to predict a geometry of higher energy because the landscape does not offer any further minimizing solutions. The final point on each CG plot merely shows this final step. The penultimate step is then treated as the final optimized solution. The 0 K plots are truncated to 500 out of the 5000 completed timesteps since the molecular energy does not fluctuate by more than 0.0001 kJ mol −1 following this point. Interestingly, the CG runs almost reach energy convergence on a similar time scale to 0 K. The lack of a clear difference between the number of timesteps is unexpected but can be explained by the fact that such a molecular system could be too small to really benefit from the CG approach. Soon-to-be published work will expand this approach to larger molecules where the number of time steps required for a successful optimization can be readdressed and confirmed with scaled-up examples.

Figure 4

T500 molecular model geometry optimization trajectory steps with SP1 (blue), SP2 (red) and SP3 (green) starting points: ( a) Set 1 (0 K and 1 fs timestep) truncated at 500 steps where the energy fluctuation is <0.0001 kJ mol −1 and ( b) Set 3 (GC and 1 fs timestep) with no truncation. The x-axis marks the timestep number. In the left panels, the y-axes denote molecular energy; in the right panels the y-axes denote Δ E (current energy - previous energy). All energies are in kJ mol −1.

Finally, Set 1 serves as an example to discuss the geometrical aspect of the optimization results. Table 3 summarizes each model's geometries optimized using Set 1's parameters. The quantitative data are augmented by data bars representing the deviation of optimized bonds and angles with respect to the QM values. Unsurprisingly, the 100 model is the worst performing, consistently producing geometries that are the most different from the QM optimum. However, even for the 100 model, the correct bond stretches are reproduced to within 0.01 Å, whereas the H-O-H angle is less accurate. The 300, 500 and T500 models all perform relatively similarly, with bond stretches within 0.007 Å and angles within 0.39°. The TE500 model performs exceptionally well, predicting both bond stretches to within 0.0005 Å and the angles to within 0.09°. Note that the TE500 model is the only model to reproduce a symmetrical final geometry. TE500 also optimizes to almost the exact same geometry for all three starting points (with Δ E = 0.14 kJ mol −1), with the final angles differing < 0.001°. All other models return different optimized angles for each SP, with some models reporting the same bond stretches across some SPs. The lack of consistency for the bond angle confirms that angular features are less energetically influencing than radial features (bond distances).

The above analysis is presented from a critical point of view in order to properly scrutinize the results. However, as seen from the energy results in Table 2, we are often working within the accuracy threshold, and the error margins seen across all the results presented in Table 3 are very low for all but the 100 model. Thus, like the energy analysis, the geometrical features are proven to optimize to their correct values, within very small error margins.

Starting from a geometry outside the training set

Here we report on the robustness of the models when the optimization is initialized from starting point geometries (SPs) with energies that are all outside of the training range. The set of 4 starting points will be referred to as SP-OUT 1, SP-OUT 2, SP-OUT 3* and SP-OUT 4*, which start with the following Δ E values: 195.24 (+6.53), 201.15 (+13.45), 300.50 (+111.79) and 592.21 (+403.5) kJ mol −1, respectively. The geometrical features of each SP-OUT system may be found in Table 4. The asterix, *, indicates that this starting point contained geometric features outside of the training range, and the bracketed values represent the difference in energy from the maximum of the trained energy range, in this case 188.71 kJ mol −1 for the T500 model. It was important that at least some of the geometric features lay within the training set range, otherwise the model would not be expected to perform and produce any relevant results. SP-OUT 4 is the only starting point lacking any geometrical features within the training set range, but is included for comparison. To reiterate, when making predictions, a kriging model will default to the mean value ( μ) when correlation between the training features and an example point's features vanish, i.e. when a geometry drifts too far outside of the training range. In the case of water, when one (or two) geometric feature(s) are in this position, the remaining two (or one) feature(s) are responsible for guiding the molecule back to within the operable training range. Currently, it is unknown to what extent this is possible. Should all three geometric features be outside the training range, the model is expected to fail and give a poor final geometry.

Table 5 reports the optimized energies from the four SP-OUT runs, again using the same four parameter sets (Sets 1 to 4). For consistency with Fig. 4, the T500 model is used for this analysis, however, any of the other 500, T500 or TE500 models would have been suitable. Remarkably, the optimizations run successfully for each of SP-OUT1, SP-OUT2 and SP-OUT3 cases, reaching Δ E values of less than 0.14 kJ mol −1.

The final geometric data for the runs of Set 1 are given in Table 4, which is analogous to the format of Table 3. Geometrical features appear good, matching the energy optimization for all except SP-OUT4. SP-OUT 4 fails by incorrectly predicting the O1-H2 bond by +0.558 Å and the H-O-H angle by −15.138°. Elongation of the O1-H2 bond causes the O1-H3 bond to shorten and finish with a reasonable final length (only +0.008 Å from the target value). Examples of the geometric trajectory can be found in the Supplementary Information. Figures S2 and S3 depict the fluctuation of the geometric features for SP-OUT 2 Sets 1 and 4 respectively, recovering from outside of the training range and producing a good final geometry. However, Figures S4 (Set 1) and S5 (Set 4) illustrate how such a recovery never occurs with SP-OUT 4. Here, the optimizations eventually terminate with all or some of the final geometric features still not, or never, in the training set range. The behaviour of the final few steps of the optimization in Figure S5 is comparable to that observed in the failed case of [SP1/500 model/"Set 4"] (see Table 2). This behaviour characterizes an evolution outside the operable range of the models.

Scanning the landscape

As a final check of the model's robustness, a comprehensive analysis was set-up by taking, as starting points, each of the 2000 geometries generated by the distortion method. Optimization runs were carried out for each geometry with the 0 K method (to ensure consistency and proper comparison) for 2000 timesteps of 0.5 fs. From the previous investigations, these parameters were deemed suitable to obtain a bird's eye view of the set's general behaviour and detect outliers. Again, the T500 model was selected for this analysis. The energy evolution of every trajectory was then extracted and aggregated, to be plotted in Fig. 5 as differences with respect to the QM minimum. By the 2000 th step, the average energy difference reached was 0.056 kJ mol −1 (standard deviation: 0.046 kJ mol −1, minimum and maximum: 0.003 and 0.24 kJ mol −1, respectively).

Figure 5

Performance of the T500 model using the 0 K optimization: ( a) aggregated plot of the molecular energy evolution in time for each of the 2000 starting geometries considered (runs are coloured from dark to light blue allowing tracking); ( b) magnified energy evolution between the 900 th and 1000 th timesteps; ( c) distribution of energies at the 1000 th timestep, relative to the ab initio energy.

The generality of the behaviour described earlier for the three individual starting points can also be observed in the top left panel of Fig. 5: every trajectory's energy monotonically drops within ~15 kJ mol −1 of the QM minimum in less than 100 timesteps. The more slowly converging trajectories seem to mostly originate in low-lying energy starting points (darker blue). The final steps, as seen in the top right panel of Fig. 5, reveal several trajectory bundles approaching the expected QM energy minimum by a different amount (again, without clear separation of the trajectories originating from low and high-lying energy starting points). Within the 2000 optimization steps, more than 25% of the set converges within 0.01 kJ mol −1, while a bigger portion of the set (~50%) clusters around convergence within 0.05 kJ mol −1. Finally, 100% of the set converges within 0.25 kJ mol −1. The most probable reason for this behaviour is the inevitable presence of noise (at least with the current method) in the kriging model, caused by the underlying accuracy of the IQA calculations, and leading to small spurious local minima around the global minimum and a less smooth PES surface. In any case, the accuracy threshold we operate within is reasonable enough to consider the whole of the sample set reasonably converging. In order to provide a fuller picture of the convergence behaviour, adequate tools for a proper analysis of the geometry evolution are in development and results will be featured in forthcoming publications featuring a larger variety of systems.

Individual preferences of the , , and atomic energies

Here we report on the individual tendencies of each of the three IQA energies (, and ), which are used in four of the five molecular models tested. Observing the 'ideal' behaviour of each gives us an insight into the interplay between the knowledgeable topological atoms occurring in each of the optimization runs. Here, we systematically eliminate, in turn, each of the other two energy contributions. First, the and models are switched off to observe only the behaviour of in an optimization. Likewise, the next run switches off the and components to obtain the behaviour of just . The final run completes the analysis by switching off and to observe only the behaviour.

Once more, the combination of the T500 model and Set 1 was selected for these examples. At the QM minimum, all three IQA components are optimally balanced: by optimising from the corresponding QM minimum geometry we can then observe each of the individual IQA components' preferential drift. Figure 6 illustrates the resulting geometries, alongside the QM initialization geometry for reference. Table 6 accompanies the results, analyzing the geometries from Fig. 6 quantitatively, similar to Table 4.

Figure 6

Single-energy optimized water geometries using the individual , and energies. Initialization geometry is the QM minimum, and the optimizations are performed using the T500 model with parameter Set 1.

Table 6: Geometrical data for the single-energy optimized runs, using Set 1 associated to the T500 model. Optimized values are reported as relative to the QM, i.e. the value of [Resulting Feature - QM], and plotted as a relative data bar. Optimized values are reported as relative to the QM, i.e. bond distances and angles are plotted as "relative data" bars where red indicates a lower value, blue a higher value. The magnitude of each bar is marked by its length, normalized using all resulting bond distances across all three SPs. The largest bar (red, V x AA', O1-H2) is set to one unit of length. The angles are treated similarly, with the unit length bar being "blue, V cl AA'".

First, the Coulomb-only optimization () slightly compresses the O-H bond lengths but maximizes the H-O-H angle to 172.6°, indicative of the, admittedly very low, electrostatic repulsion between the two hydrogen atoms. Second, the intra-atomic-only optimization () significantly elongates the O-H bonds and makes the H-O-H angle much more acute (~50°). Unpublished results have shown that can generally be correlated with atomic volumes: where becomes lower (more stabilized) atomic volume increases and vice versa. Hence, elongation of the O-H bond lengths reflects an increase in both the oxygen and hydrogen atomic volumes. However, the H-O-H angle becoming more acute tells us that maximizing the oxygen's volume must more than compensate for the additional destabilization resulting from the mutual compression of the two hydrogens. Third and finally, the exchange-only optimization () expectedly compresses the O-H bond lengths upon optimization. This compression corresponds to a lowering of the interatomic exchange energy. This energy is numerically dominated by the O-H exchange energy. Thus, the O-H compression corresponds to a strengthening in the covalent bonding in water. However, at some point in the compression the exchange stabilization reaches an extremum, and the bond becomes increasingly unstable again. We can also deduce that the through-space interaction between the H...H atoms must not be significant since it is barely able to make the H-O-H angle more acute (~90°).

The geometries presented in this section should not be overinterpreted, as they are driven outside of the training range (in "no-man's-land") and are likely to converge due to the loss of kriging correlation: these analyses are relevant for the geometric drift tendencies only.

The analysis in this section reminds us how there is not only a complex interplay between the type of atomic energies used within the optimization, but also a balance being reached between the atoms when only a single type of IQA energy is used. The interpretation of significance48 of each IQA energy could benefit from the perspective provided by such optimizations, at least at the level of chemical intuition, as it gives some insight about where each energy component "pushes" the molecule to go toward. No such analysis is possible for the TE500 model, for which the earlier results (see Section 3.2) already illustrate the preferential behaviour of the energy for an atom in this water model.

These results stand in contrast to those that would be made on a traditional force field potential. For a single water molecule, a typical bonded potential would consist of a single harmonic angle-bend and a pair of identically-typed harmonic bond-stretch terms. Turning off one of these terms while maintaining the remaining two will result in a completely unphysical potential. Removing the angle-bend term gives a potential where the H atoms can occupy the same position with no energetic penalty, while removing either bond-stretch term allows the corresponding O-H interatomic distance to take any value. Neither of these altered potentials can provide any physical insight, in stark contrast to the potential described herein, emphasizing the lack of correspondence between standard force field bonded energy terms and their underlying quantum origins.

Discussion

Contextualization of FFLUX

A thorough dissection of the force field research literature proves that the FFLUX methodology is unique and indeed novel. This approach goes further than multipolar force fields (such as SIBFA and AMOEBA), which in turn innovate the traditional point charge force fields (such as AMBER).

The work we presented here is a very detailed theory-versus-theory assessment rather than a theory-versus-experiment comparison. This means that we set up a novel computational scheme (FFLUX) that is asked to reproduce, as best as it can, the original quantum mechanical data that it was trained for. We carefully demonstrated a proof-of-principle, at some level of theory that is not the best possible because it does not have to be the best possible, that is, for the current purpose of proof-of-principle.

We worked with Hartree-Fock because the IQA energy contributions are very clear in this case. We took a small basis set because it compensates the inherent error introduced by the limitations of the Hartree-Fock Ansatz. As long as the proof-of-principle of the machine-taught topological atoms reaching the global energy minimum together is solid we have reached our goal, no matter at which level of theory this is achieved. The current article reports the first ever such study, and this successful proof-of-concept opens an exciting avenue for a host of applications, on larger molecules and complexes, eventually reaching condensed matter simulations via a rigorous bottom-up research program.

Comparing to other traditional force fields such as AMBER or CHARMM is also moot at this stage. The geometry of a single water molecule in a force field method can be optimized by simply looking up the values for the two reference bond lengths and the one reference angle in the appropriate parameter set. Note that there will be no non-bonded interactions in a standard force field for a single water. These values are almost always chosen by comparison to quantum chemical results anyway; how well a traditional force field describes a single water molecule is of very limited interest. Because of the novelty of the FFLUX method there is already enough to explain on a single water molecule. This system illustrates in an "uncluttered" way how much parameterization effort is already involved. This article raises points that are salient to proving that the method can be used in much more complicated applications (which are the subject of future publications).

We do not intend to ignore the vital question of how a "force field makes contact with experiment". Of course the long-term goal is to make reliable predictions that are experimentally verifiable. More excitingly, we eventually intend to use FFLUX in the area of nucleation, which experiment cannot probe in its very early stage. This is an example where FFLUX's reliability will be crucial.

In terms of comparison with experiment at this current early stage, the relevant experimental data for a single water molecule are then the 0 K geometry (and potentially its vibrational frequencies), since the geometry is what we aim to compute. However, this geometry is actually not known experimentally: bond lengths are never known as (quantum mechanical) r e values but as any one of derived (i.e. treated) bond lengths (e.g. r z or r g or r α ...). In our work we use the Hartree-Fock method as a surrogate for experiment, being straightforwardly (if expensively) replaceable with any other quantum mechanical method in the construction of FFLUX models. This successful and safe strategy defers the question of comparison to experiment to later applications where, for example, liquid density can be evaluated from simulation and directly compared to experiment.

Our overall strategy has always been bottom-up: start from small systems (such as a single water molecule) and then upscale. In the current proof-of-concept stage it is important and indeed sufficient to compare FFLUX with first-principles (i.e. quantum mechanical) data, which is what is done in this paper. In later work we will carry out simulations on liquid water (as we have done with rigid body water molecules, equipped with multipolar electrostatics). We will then, as we have done in our previous publications, compare the computed structure and dynamics with neutron diffraction data, and kinetic and thermodynamic quantities (e.g. self-diffusion coefficient, isothermal expansion coefficient, heat capacity at constant pressure etc.) with measured values.

Compared to a traditional force field, a FFLUX model will always take more computational time but then the former provides less information than FFLUX; indeed, one should not compare profoundly dissimilar objects. FFLUX is a force field that "sees the electrons", which traditional force fields do not. As demonstrated in the last part of the article, FFLUX is aware of the internal energy of an atom (i.e. intra-atomic), the electrostatic interaction energy and the exchange energy. These are quantum mechanical data that originate from the wave function itself. Traditional force fields do not contain this information. Moreover, FFLUX also stores atomic multipole information (not active in the optimization of water because of divergence of the multipole expansion but the electrostatic interaction is covered by the IQA term V cl anyway), and also stores their polarization more importantly, which is absent in traditional force fields. Again, this extra information adds to the cost of FFLUX.

It has been clear from the start of the FFLUX project, that FFLUX will be a more expensive force field computationally but, at the beginning of the FFLUX project many years ago, we aimed for computers of the near future, which now exist. Each year passing enables more expensive calculations to become feasible, and thus the reliability and accuracy of FFLUX will increasingly benefit the systems that are currently only within the application radius of traditional force fields.

Building up a database of transferable models will be a computationally expensive task but will remain a one-off, not burdening the user who wants to be shielded from this activity. However, there is nothing stopping users adding to a database of models, much like the PDB or CSD. Building models (more specifically, obtaining the IQA-QM data) is the expensive stage. However, the high transferability of our topological atoms has already proven to reduce the necessary workload.

CPU timings for a single water optimization can be reported at this stage although the code has not been optimized yet. Without compiler optimisation and in debug mode on a single Intel® Core™ i5-2410M CPU@2.30 GHz processor, DL_POLY takes, for 2000 steps (of the iterative geometry optimisation) and 500 training examples, using all three energies, 1.92 seconds, or about one millisecond per step. This result is very encouraging because FFLUX already defeats the original ab initio calculation by at least three orders of magnitude (and this is just for Hartree-Fock, which is a cheap ab initio method), prior to any source code optimization having taken place.

Finally, we point out that in terms of accuracy the method relies on much more representative physics than a harmonic force field.

Final Considerations in connection with Future Work

Looking at the nature of the FFLUX approach and its future brings up five topics that benefit from some extra comment at this point in time.

First, many-body effects often feature in force field discussions, particularly in the treatment of liquid water in terms of (long range) perturbation theory. While focusing on water clusters, the concept of many-body effects determines to what extent a water trimer for example, can be described, energetically, by three water dimers. We appreciate this type of analysis but it does not directly affect the IQA analysis that we use in FFLUX. The key point is that all IQA quantities are always extracted from the full wave function, involving all water molecules in the cluster. In other words, the electron density is obtained through the Self-Consistent-Field (SCF) procedure, which relaxes all the orbitals with respect to one another, and thereby automatically accounts for many-body effects. We also note that FFLUX incorporates polarization through kriging models for atomic multipole moments responding to their precise environment. Unlike other polarizable force fields, FFLUX does not invoke and a SCF procedure during the molecular dynamics production run. The present method aims, first and foremost, at improving the realism of molecular modelling and bringing it closer to quantum mechanics, which is known to have an even higher computational cost. However, the model offers optimisation opportunities, e.g. in terms of mass-parallelism, with the prospect of reaching quasi-linear scaling on distributed hardware.

Secondly, there are advantages and disadvantages to FFLUX. Attractive features encompass (i) the near ab initio accuracy without having to carry out the ab initio calculation in the "production phase", (ii) much faster performance than that of an ab initio calculation, especially if the latter is carried out with a high level of theory, (iii) robust chemical insight, physically rooted and free from assumptions. For example, one can find out why an energy minimum exists in terms of the balance between intra-atomic self-energy, electrostatic and exchange(-correlation) energies, (iv) no need for a penetration correction, (v) no need for damping functions (because there is no polarisation catastrophe), (vi) FFLUX is well grounded in the literature with precise answers regarding its various aspects, (vii) the topological partitioning is parameter-free and reference-free, (viii) diffuse Gaussian functions are not problematic, and finally (ix) Kriging handles high-dimensional feature spaces well with a relatively small number of data points. Amongst the current disadvantages one can think of (i) the computational expense of IQA, (ii) the performance of kriging models depending on atomic integration errors, (iii) an improved sample point selection approach to correctly represent the configuration space at hand, and finally (iv) the challenge of defining ALFs in condensed matter, in the presence of significant intramolecular hydrogen bonds or π-π stacking, and when molecules travel over large distances.

Thirdly, in terms of suitability for quantum dynamics simulation, FFLUX can in principle be applied to adiabatic Born-Oppenheimer ab initio MD, as it provides the ground state total energy and forces for a chemical system required to evaluate the forces on the nuclei. Whether or not this is to be viewed as quantum dynamics is open to interpretation due to the use of the machine learning method as an intermediary between nuclear coordinates and the total energy. Whilst definitely not a force field method, the described potential is in some sense empirical. Beyond ab initio MD, FFLUX is not useful in its current state for non-adiabatic calculations as the system wavefunction is only implicit in the output of the machine learning models. Molecular orbital or electron density information is not produced by the Kriging models; the method shares an inherent ground-state nature with force fields. However, for reactions or systems with small band gaps, one would certainly need to be able to incorporate non-adiabatic effects.

Fourthly, FFLUX can be applied to bond breaking, in principle, including chemical reactions and changes in metal coordination. A kriging model can be taught potentially any behaviour, given the correct data. Note that in this study we use kriging to find relationships between energies and geometries, with the geometries being represented by features. Those features can take any value and be defined as arbitrarily as the user requires. Thus it is not only possible to tackle a problem such as bond breaking, it should require no additional terms or fixes such as those found in several force fields when attempting to model new phenomena.

Fifthly and finally, the number of features in the current work is very small but kriging scales exceptionally well. In the recent past we have tackled over 100 features for larger molecules18,49 and molecular clusters while maintaining very good kriging predictions.

Conclusion

For the first time, atomic kriging models have been "set in motion", through the associated (analytical) forces. Geometry optimizations have been successfully carried out, yielding energies and geometries in agreement with the QM optimum. While used as the model's seed, the latter is not part of the training set: the ability of our kriging atomic models to generally reproduce a sampled molecular potential is then fully confirmed.

A variety of kriging models were analyzed and compared: complementing the picture provided by S-curves, optimization stands as a new validation tool, closer to practical purposes and sensitive to gradient prediction errors. As expected, models including more training points yield optimized structures closer to the QM reference, both in energy and geometry. The robustness of the kriging models is demonstrated by their ability to fall back into their optimum even when starting outside of their conformational training range. Initializing the optimization from every generated sample point did not reveal any major spurious minimum or other shortcoming on the potential energy surface, thereby corroborating the kriging method.

Finally, chemical insight provided by the IQA energy decomposition is preserved through our method, where observations consistent with intuition have been made by isolating individual contributions.

Further encouraging results are soon-to-be published, featuring more complex molecules and a deeper analysis of the potential energy surfaces, in particular the forces and the agreement between energy predictions and original QM energy.

Additional information

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. 10.
Krige, D. G. Two-dimensional weigthed moving average trend surfaces for one valuation. J. South African Inst. Mining Metall. 13-38 (1966).

2. 15.
Davie, S. J., Di Pasquale, N. & Popelier, P. L. A. Incorporation of Local Structure into Kriging Models for the Prediction of Atomistic Properties in the Water Decamer. J. Comput. Chem. 37, 2409-2422 (2016).

3. 25.
Maxwell, P. I. & Popelier, P. L. A. Accurate prediction of the energetics of weakly bound complexes using the machine learning method kriging. Struct. Chem. in press doi: https://doi.org/10.1007/s11224-017-0928-9 (2017).

4. 26.
Cardamone, S. & Popelier, P. L. A. Prediction of Conformationally Dependent Atomic Multipole Moments in Carbohydrates. J. Comput. Chem. 36, 2361-2373 (2015).

5. 29.
Frisch, M. J. et al. GAUSSIAN09, revision B. 01. Gaussian Inc., Wallingford, CT, USA (2010).

6. 30.
Malcolm, N. O. J. & Popelier, P. L. A. The full topology of the Laplacian of the electron density: scrutinising a physical basis for the VSEPR model. Faraday Discussions 124, 353-363 (2003).

7. 37.
AIMAll, T. A. K. TK Gristmill Software, Overland Park KS, USA, (aim.tkgristmill.com) (2016).

Download references

Acknowledgements

P.L.A.P. acknowledges the EPSRC for funding through the award of an Established Career Fellowship (grant EP/K005472).

Author information

Present address: Joint BioEnergy Institute, Emeryville, CA, 94608, USA

Affiliations

1. Manchester Institute of Biotechnology (MIB), 131 Princess Street, Manchester, M1 7DN, United Kingdom

François Zielinski

, Peter I. Maxwell

, Timothy L. Fletcher

, Stuart J. Davie

, Nicodemo Di Pasquale

, Salvatore Cardamone

& Paul L. A. Popelier

2. School of Chemistry, University of Manchester, Oxford Road, Manchester, M13 9PL, United Kingdom

François Zielinski

, Peter I. Maxwell

, Timothy L. Fletcher

, Stuart J. Davie

, Nicodemo Di Pasquale

, Salvatore Cardamone

& Paul L. A. Popelier

3. Biomass Science and Conversion Technology Department, Sandia National Laboratories, Livermore, CA, USA

Contributions

F.Z. wrote the new subroutines in DL_POLY and conducted all tests assisted by P.I.M., T.L.F. provided the GAIA code, S.J.D. provided overall consulting, N.D.P. wrote the kriging code and S.C. the distortion code, and M.J.L.M. provided the implementation of analytical forces. P.L.A.P. conceived, guided and micromanaged the whole project, and also rewrote large parts of the original manuscript.

Competing Interests

The authors declare that they have no competing interests.

Corresponding author

Correspondence to Paul L. A. Popelier.

Electronic supplementary material

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

# data[].attributes.attachment_link.title
Geometry Optimization with Machine Trained Topological Atoms

# data[].attributes.attachment_link.description
Altmetric: 1 Citations: 1 More detail Article | Open Scientific Reports 7, Article number: 12817 (2017) doi:10.1038/s41598-017-12600-3 Download Citation Received: Accepted: Published online: The geometry optimization of a water molecule with a novel type of energy function called FFLUX is presented, which bypasses the traditional bonded potentials.

# data[].attributes.headline
MOHAMAD SHAFIQ BIN MOHD SUKRI 181872              

# data[].attributes.subject
MOHAMAD SHAFIQ BIN MOHD SUKRI 181872              

# data[].attributes.body
Computational modeling of drug transport across the in vitro cornea

Pak, J., Chen, Z. J., Sun, K., Przekwas, A., Walenga, R., & Fan, J. (2018).

- The cornea presents a unique barrier design and opportunity for investigation using Q3D modeling

- A novel quasi-3D (Q3D) modeling approach was developed to model an interconnected stack of different corneal tissue layers with membrane junction conditions assigned between the tissues.

- Lipophilic Rhodamine B and hydrophilic fluorescein were used as drug analogs.

- The model predictions for both hydrophilic and lipophilic tracers were able to match the experimental measurements along with the sharp discontinuities at the epithelium-stroma and stroma-endothelium interfaces.

# data[].attributes.attachment_link.title
Computational modeling of drug transport across the in vitro cornea - ScienceDirect

# data[].attributes.attachment_link.description
A novel quasi-3D (Q3D) modeling approach was developed to model networks of one dimensional structures like tubes and vessels common in human anatomy such as vascular and lymphatic systems, neural networks, and respiratory airways. Instead of a branching network of the same tissue type, this approach was extended to model an interconnected stack of different corneal tissue layers with membrane junction conditions assigned between the tissues.

# data[].attributes.headline
KU AIDA SUZIANA BT KU JAMALUDDIN (182227)

# data[].attributes.subject
KU AIDA SUZIANA BT KU JAMALUDDIN (182227)

# data[].attributes.body
Title: The impact of protein interaction networks' characteristics on computational complex detection methods.
Author: Xiaoxia Liu, Zhihao Yang, Ziwei Zhou, Yuanyuan Sun, Hongfei Lin, Jian Wang and Bo Xu.
Year: 2018
Date available online: 5 December 2017

What the research is about:

*Analyze the biological and topological characteristics of multiple PPI networks.
*Compare the performances of six classic computational complex detection methods on the protein interaction networks.
*Relations between the complex detection method and characteristics of PPI datasets (which may be useful to improve the performance in protein complex detection area).
*Change topological characteristics of PINs by adding protein characteristics of PINs with performances of protein detection methods.
*Compare the complexes detected by two different methods and the real complex in order to clearly describe the difference between two methods as well as locate the proteins not detected by either method.

Computational Method use:
1) MCODE (Molecular Complex Detection)
2) MCL (Markov Clustering): to detect functional modules and protein complexes.
3) CMC (Clustering based on merging Maximal Cliques)
4) RRW: discover protein complexes within large-scale protein networks.
5) ClusterONE (Clustering with Overlapping Neighborhood Expansion): predicts overlapping protein complexes.
6) COACH: detects protein complexes based on core-attachment structure of complex in two stages.

URL:
https://www.sciencedirect.com/science/article/pii/S0022519317305362

PHOTOS:

# data[].attributes.headline
The use of quantum chemical methods in corrision inhibitor studies

# data[].attributes.subject
The use of quantum chemical methods in corrision inhibitor studies

# data[].attributes.body
Gokhan Gece
27 August 2008

Through this method, it can be used in searching for compounds with desired properties using chemical intuition and experience into a mathematically qualified and computarized form. A set of mathematical equation are capable of representing accurately the chemical phenomenon when there is a correlation between the structure and activity of property is found.

URL:

Naemah Solihah Bt Jonaidi
182243

# data[].attributes.attachment_link.title
The use of quantum chemical methods in corrosion inhibitor studies - ScienceDirect

# data[].attributes.attachment_link.description
Quantum chemical methods are particularly significant in the study of electrochemistry and provide researchers with a relatively quick way of studying the structure and behaviour of corrosion inhibitors. The originality of this review article is based on the fact that it is the first and unique general reference for all those interested in the use of quantum chemical methods in corrosion inhibitor studies.

# data[].attributes.body
Electron Affinity of Liquid Water
Alex P. Gaiduk, Tuan Anh Pham, Marco Govoni, Francesco Paesani & Giulia Galli
08 September 2017
doi: 10.1038/s41467-017-02673-z

- Prediction of electron affinity of liquid water for better understanding about the mechanisms of redox reactions in aqueous system.
- Methods are the first principles which combined path-integral MD simulations with ab initio potentials, and state-of-the-art electronic structure methods based on MBPT.

URL :

# data[].attributes.attachment_link.readable_text
Altmetric: 102

More detail

Article | Open

Nature Communications 9, Article number: 247 (2018)

doi:10.1038/s41467-017-02673-z

Download Citation

Received:

Accepted:

Published online:

Abstract

Understanding redox and photochemical reactions in aqueous environments requires a precise knowledge of the ionization potential and electron affinity of liquid water. The former has been measured, but not the latter. We predict the electron affinity of liquid water and of its surface from first principles, coupling path-integral molecular dynamics with ab initio potentials, and many-body perturbation theory. Our results for the surface (0.8 eV) agree well with recent pump-probe spectroscopy measurements on amorphous ice. Those for the bulk (0.1-0.3 eV) differ from several estimates adopted in the literature, which we critically revisit. We show that the ionization potential of the bulk and surface are almost identical; instead their electron affinities differ substantially, with the conduction band edge of the surface much deeper in energy than that of the bulk. We also discuss the significant impact of nuclear quantum effects on the fundamental gap and band edges of the liquid.

Introduction

Upon exposure to light or ionizing radiation, liquid water may undergo a series of transformations involving electronic excitations, ionization of solvated species, and formation of free radicals and solvated electrons1,2,3,4,5,6,7. These processes occur, for example, in (photo-)electrochemical cells8,9, biological molecules10, as well as in atmospheric water11. Many of them are governed, at least in part, by the nature of the electron transfer and binding energies in aqueous solutions. Hence a detailed understanding of the interaction of light with water, and in general of redox reactions in aqueous systems, requires a precise knowledge of the liquid ionization potential (IP) and electron affinity (EA).

Direct measurements of the IP of liquid water are available from photoelectron spectroscopy, obtained by linear extrapolation to zero energy of the lowest binding-energy peak (1 b ₁). The measurement of Delahay et al.12 yielded 10.06 eV for the water IP, while more recent photoelectron spectroscopy studies led to a slightly smaller value of 9.9 eV13, leading to an overall agreement that the water IP is ~10 eV. On the other hand, direct measurements of the liquid EA, or the energy gain due to the injection of an electron into the liquid, are not available, due to the rapid solvation of electrons in water1,2,14. Current estimates of the EA are based on thermodynamic arguments involving adiabatic and vertical ionization energies of the solvated electron7,15,16, computed reorganization energies, and specific assumptions on the behavior of photoionized electrons in water. As a result, the reported values of EA span a wide range between 0 and 1.2 eV.

The only direct measurement of the EA of condensed aqueous systems was performed for the surface of a thin film of amorphous ice (a-ice) by pump-probe femtosecond spectroscopy20,21. Ice films only 2.5-bilayer thick were grown on a copper electrode, hence the measurement probed the energy of the lowest empty state of the aqueous surface. This study revealed a short-lived (~40 fs) spectral feature at −0.8 eV below the vacuum level, which was attributed to an electron residing at the bottom of the conduction band of amorphous ice interfaced with copper20.

Multiple predictions of water IP and EA from quantum calculations produced conflicting results depending on the structural models and electronic-structure methods employed. For example, many-body perturbation theory (MBPT) calculations using wavefunctions derived from density functional theory (DFT) and the PBE approximation22 yielded a value of 0.7 eV for the EA (configurations were generated either with classical potentials or the PBE functional); this value appeared to be consistent with estimates reported in the literature15. The quasiparticle self-consistent GW calculations of Chen et al.25, who generated a water structural model using van der Waals functionals and included nuclear quantum effects, obtained EA = 0.5 eV. A later thermodynamic integration study performed by the same group7 using hybrid functionals23 without nuclear quantum effects, predicted a larger value of 0.97 eV.

Here we present accurate predictions of the EA of liquid water and its surface using a combination of path-integral molecular dynamics with ab initio potentials, and MBPT calculations, thus determining a key missing property of electrons in water. We compared our results for the water surface with experiments on a-ice, thus validating our computational strategy against the only available direct measurement. We then used the same protocol to calculate the EA of the bulk, which we predict to be 0.1-0.3 eV. This value is similar to that estimated by Coe et al.16, although we show below that the agreement is fortuitous and effectively derives from inaccurate assumptions made in ref.16. We also revisit several other experiments and the assumptions that led to current literature values of the EA.

Results

Computational model

The structural models of liquid water were generated using path-integral molecular dynamics (MD) simulations with interatomic interactions described with the many-body MB-pol potential energy function27,28,29, which has been shown to accurately predict the properties of water in gaseous and condensed phases30. Electronic-structure properties were then computed for MD snapshots utilizing MBPT at the G ₀W ₀ level of theory, starting from eigenvalues and wavefunctions determined from hybrid density functional calculations. The water IP was computed by extrapolation of the (1 b ₁) peak energy, similar to the procedure followed to analyze experimental data13; the EA was computed as the lowest unoccupied single-particle energy level of the liquid (see SI). Within MBPT, the energies −IP and −EA correspond, respectively, to the valence band maximum (VBM) and conduction band minimum (CBM), which were referred to vacuum using a procedure commonly utilized to compute band offsets in semiconductors (see Supplementary Methods).

In order to validate our computational protocol, we first compared the energies computed for the occupied electronic states of water with well-established liquid jet experimental data13,36. Consistent with our previous studies37,38, we found that the most accurate protocol for computing the binding energies of electrons in water is the G ₀W ₀ approximation based on dielectric-dependent hybrid (DDH) functionals39,40. In particular, G ₀W ₀ calculations with the range-separated hybrid (RSH) yielded a mean absolute error (MAE) with respect to measured spectra of only 0.18 eV for trajectories including nuclear quantum effects (NQE) and 0.38 eV for classical trajectories.

In addition, to assess the performance of the same level of theory for the unoccupied states of the liquid, we relied on the comparison of our results with highly-accurate quantum chemistry calculations of the EA of water clusters41, due to the lack of experimental data for the bulk system. Specifically, we compared G ₀W ₀ electron affinities with those obtained with the coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] method for a series of water hexamers, and found that overall, G ₀W ₀/RSH quasiparticle energies agree with the CCSD(T) results within 0.05 eV. This accuracy is comparable to the statistical uncertainty of our calculations.

Energy levels of bulk water and the water surface

Figure 1 (left-hand side) summarizes the computed VBM and CBM of a slab representing the water/vacuum interface, as obtained at the G ₀W ₀/RSH level. The computed CBM values of −0.88 eV (G ₀W ₀/RSH) and −0.79 eV (G ₀W ₀/sc-hybrid) are in agreement with the direct measurement of Stähler et al.20, who reported a conduction band minimum of −0.8 eV at the interface of amorphous ice and a metal electrode. The agreement with experiment found here for the a-ice surface, expected to be representative of that of the liquid as well, further proves that our computational protocol accurately predicts the energy of the empty states of water.

Fig. 1

Computed electronic energy levels of liquid water. Upper panel: Positions of the valence band maxima (VBM, blue) and conduction band minima (CBM, red) of the surface of the water and bulk water computed using the classical and path integral molecular dynamics with the MB-pol potential. All values are in eV. The energy levels were computed using G ₀W ₀ starting from hybrid DFT; the range given above (thick bars) corresponds to results obtained with range-separated (RSH) and self-consistent hybrids (see text and Supplementary Tables 7- 10). Lower panel: Snapshots of the surface of water (left) and of bulk liquid water (right) from path integral MD simulations, together with isoprobability contours (set at 40%) of the lowest unoccupied molecular orbital, as obtained using the RSH functional. The interface between water and vacuum (gray isosurface in the lower left panel) is pictorially represented using a constant density surface defined following the procedure reported in ref.69

We now turn to discussing the EA of the bulk liquid. As Fig. 1 shows, the position of the VBM of the surface and bulk water are similar (differing by less than 0.1 eV); however the position of the bulk CBM is ~0.6 eV higher than that of the surface. Our results obtained with G ₀W ₀/RSH (−0.29 eV) and G ₀W ₀/sc-hybrid (−0.17 eV) yield EA of the bulk within 0.1 and 0.3 eV, when including statistical uncertainties of our calculations (0.04-0.05 eV). This range has been obtained using MD trajectories inclusive of nuclear quantum effects; with classical trajectories and the same level of electronic structure theory, we found smaller values: 0 to 0.1 eV.

Nuclear quantum effects

The results presented in Fig. 1 highlight the influence of nuclear quantum effects on the electronic structure of liquid water. In particular, the inclusion of NQE affects the VBM and CBM position by 0.3 and 0.2 eV, respectively (the energy of the VBM is increased and that of the CBM lowered) yielding a gap 0.5 eV smaller than that calculated for classical trajectories. The decrease in the gap is consistent with the presence of longer O-H bonds (0.99 ± 0.14 Å) when NQE are included (to be compared with 0.97 ± 0.05 Å for classical trajectories): Indeed, ref.22 showed that water models with longer O-H bonds exhibited smaller gaps. Accounting for NQE also results in a wider bond-pair peak in the distribution of maximally localized Wannier centers (Supplementary Fig. 6), consistent with larger O-H bond-length fluctuations.

Our result for the fundamental gap is in agreement with the study of Chen et al.25 reporting a 0.5-0.7 eV reduction of the band gap of water due to nuclear quantum effects. Instead our findings are inconsistent with the density of virtual states reported in ref.34 showing a rather large shift of the CB edge upon inclusion of NQE (of almost 1 eV). We also note that the fluctuations of the VBM and CBM observed in our quantum trajectories are of the order of 0.4 and 0.1 eV for the bulk and 0.6 and 0.4 eV for the surface, indicating substantial fluctuations of band edges depending on the nuclear configuration, which should be taken into account when modeling chemical reactions at aqueous interfaces.

The marked difference in the EA of the bulk and the surface is reflected in the localization properties of their respective conduction band edges. In Fig. 1 we report the isosurface of the wavefunctions corresponding to the CBM of the surface and the bulk, showing the different localization properties of the two states. Consistent with previous reports42,43, we found that the bulk unoccupied edge is delocalized over the entire supercell, while that of the slab is localized at the surface, in proximity of broken hydrogen bonds. The different bonding environment of the surface, with respect to that of the bulk, is reflected in the values of effective polarizabilities44: As expected from our previous work on salts45, we found that the effective molecular polarizabilities in the slab, where more broken hydrogen bonds are present, are smaller than in the bulk. In addition the distributions of polarizabilities and dipole moments are broader in water with NQE than in the liquid treated classically (Supplementary Figs. 5 and 6).

Discussion

Our prediction for the liquid vertical EA (0.1-0.3 eV) differs from data often quoted as experimental results in the literature, e.g., 0.74 eV from ref.15. However, these data are not obtained from direct measurements; rather, they are derived using both experimental and theoretical results, under specific assumptions which, as we show below, need to be revisited. By combining our theoretical prediction of the vertical EA with the most recent measurements on solvated electrons, we discuss below the construction of an energy diagram for an electron in liquid water.

The vertical (EA) and adiabatic (AEA) electron affinities of water are defined as:

and

where E water is the energy of pristine, neutral water and E that of the liquid with a non-solvated electron; E aq is the energy of the liquid with a solvated electron e − (aq). The vertical detachment energy (VDE), or energy necessary to remove an e − (aq) from water, is defined as

where E water cavity refers to the energy of the neutral liquid with a defect (cavity) created by a solvated electron. The energy required to create the defect is the reorganization energy λ:

Finally, E e - aq * is the energy of the solvated electron's first excited state with excitation energy μ =  E e - aq *  −  E aq .

Figure 2 shows the electron energy diagram in liquid water constructed from our EA calculations and the most recent measurements for AEA, VDE, and μ. Specifically, AEA = 1.34 eV was obtained by extrapolation of water cluster data46; VDE = 3.7 eV was measured by photoemission spectroscopy including corrections for surface scattering effects47; and μ = 1.73 eV is a well-established position of the maximum in the measured optical absorption spectrum of the solvated electron48. The reorganization energy λ derived from the data of Fig. 2 is 2.36 eV, larger than the most recent computational estimates of 1.99-2.17 eV7.

Fig. 2

Energy diagram of an electron in water. The energies reported in this diagram are from the most recent experimental results available in the literature and the value of the vertical electron affinity (EA) is determined in this work (0.1-0.3 eV). The vertical detachment energy of the solvated electron is VDE = 3.7 eV47 and the adiabatic electron affinity of water is AEA = 1.34 eV46. The difference between VDE and AEA is the water reorganization energy upon solvation of an electron: λ = 2.36 eV; note the difference of more than 1 eV with the values used in the left panel of Fig. 3. E e - aq * corresponds to the excited state of the solvated electron e − (aq), with excitation energy μ = 1.73 eV48. All values are in eV

The diagram obtained here differs substantially from those obtained using the values reported by Bernas et al.15 and by Coe et al.16, shown on the left-hand and right-hand sides of Fig. 3, respectively. Focusing first on the results of ref.15, the authors considered AEA = 1.47 eV, λ = 1 eV49,50, and μ = 1.73 eV48. This led to a VDE of 2.47 eV, which is significantly different from that reported by recent measurements (3.3-3.8 eV). This discrepancy originates from the incorrect magnitude of the reorganization energy used in ref.15, a value obtained from classical MD simulations49,50, which turned out to be underestimated by more than 1 eV with respect to all best estimates known today. In addition, a crucial assumption of ref.15 was that the first excited state of the solvated electron coincides with the bottom of the conduction band of water, hence E  ≡  E e - aq * . This assumption was later disproved and is now known to be incorrect54.

Fig. 3

Energy diagrams of solvated electron inferred from past work. The diagrams are derived from the work of Bernas et al.15 (left) and Coe et al.16 (right). The notation is the same as in Fig. 2. The zero of both diagrams corresponds to the energy of pristine, neutral water, E water. Bernas et al. used AEA = 1.47 eV70, while Coe et al. adopted AEA = 1.72 eV. On the left, λ = 1 eV49, 50 was taken from a theoretical study and μ = 1.73 eV from optical absorption experiments48. On the right, VDE = 3.32 eV55 was extrapolated from cluster data. Ref.15 (left) assumed that the energy of pristine water with an added electron (prior to ionic relaxation) coincides with the first excited state of the solvated electron E e - aq *

Finally, we note that Bernas et al.15 defined the bottom of the conduction band of water from the energy difference V 0 = E e-- E water cavity , which, however, is not equal to −EA. Using the energy diagram of Fig. 3, one would obtain a negative value for the EA of liquid water: −0.26 eV.

Coe et al.16 (right-hand-side diagram) used the values VDE = 3.32 eV and AEA = 1.72 eV, both derived from extrapolated experimental data for water clusters55, and obtained a larger reorganization energy ( λ = 1.60 eV) than Bernas et al.15 They then computed an upper bound to the EA as (AEA −  λ) by assuming that the energy difference ( E  −  E aq ) is always smaller than the reorganization energy. However, the derivation incorrectly assumed that the creation of a cavity in pure water implies an energy gain, rather than a cost, as pointed out in ref.7. If energy gains and costs are correctly taken into account, the only conclusion that can be drawn from the data of ref.16 is EA ≤ 1.72 eV.

We compare our results with other computational studies, focusing on first-principles calculations. Using several water models, previous studies predicted the EA of the bulk between 0.7 and 0.9 eV at the G ₀W ₀/PBE level of theory. As shown in Supplementary Table 7, we obtained a similar result (0.99 eV) for classical MB-pol water with the same electronic structure methods. Our best prediction for the bulk water EA is smaller than that reported recently using first-principles MD trajectories with van der Waals functionals and nuclear quantum effects (0.5 eV)25. These calculations used self-consistent GW with approximate vertex corrections and smaller cells (32 water molecules) than those adopted in our study (64-256 molecules).

Knowledge of the water EA is key to understanding mechanisms of redox reactions in aqueous systems that involve either molecular species or solid surfaces53,56,57, yet no direct experimental measurement is available. In this paper we reported predictions of the EA of bulk water and the surface of water, obtained entirely from first principles, by combining path-integral MD simulations with ab initio potentials, and state-of-the-art electronic structure methods based on MBPT. The accuracy of the theoretical methods used in this work was carefully checked against quantum chemistry calculations for virtual states of aqueous systems41, and against well-established experiments for electronic binding energies of the liquid. Our results for the surface were found to be in good agreement with recent pump-probe spectroscopic measurements20, thus further validating the accuracy of our computational framework.

We determined the EA of the liquid to be between 0.1 and 0.3 eV, much smaller than some of the accepted estimates present in the literature, but consistent with the speculations of Coe et al.16,18 and the measurements of the electron ejection lengths in a two-photon ionization process by the Bradforth group58,59 (The authors of refs.58,59 determined that the average ejection length of the electrons in a two-photon ionization process is roughly constant below 9.5 eV but increases rapidly after 9.5-9.8 eV, and associated this increase with the electrons being able to access the conduction band upon excitation. Given that the ionization threshold of water is 9.9 eV13, the energy of 9.5-9.8 eV corresponds to an energy level positioned ~0.1-0.4 eV below vacuum, although refs.58,59 do not explicitly mention these values for the conduction band minimum). We used our theoretical prediction of the vertical EA together with the most recent measurements for the solvated electron, to construct an energy level diagram for an electron in water, which we believe is the most accurate known to date. The diagram proposed here differs substantially from those generated using estimates reported in the literature, for which we provided a detailed, critical reassessment. We suggest that the EA of water may be determined experimentally by extending the pump-probe spectroscopic study of ref.20 to thicker a-ice films and extrapolating the results as a function of size, building, e.g., on work reported by King et al.21 on trapped electrons at interfaces.

Our results showed that while the IP of bulk water and the surface are very similar, their EAs differ by more than 0.6 eV, with the conduction band edge of the surface being much deeper in energy than that of the bulk. We also found that the band gap of water as well as both band edges are substantially affected by proton quantum effects and couplings with nuclear configurations, highlighting the importance of including energy level fluctuations when modeling chemical reactions in water. Work is in progress to investigate the same effects in aqueous solutions.

Methods

Molecular dynamics simulations

Classical and path-integral (PI) molecular dynamics (MD) simulations were performed using the ab initio MB-pol potential and a modified version of the DL_POLY code. Bulk water was simulated in the constant volume-constant temperature (NVT) ensemble at 298 K, using a Nosé-Hoover chain thermostat60 with the length of 4, and supercells containing 64 molecules. The equilibrium density ( ρ) of the liquid was determined with constant pressure and temperature (NPT) simulations, corresponding to cell sizes of 12.41 Å ³ ( ρ = 1.001 g cm −3) and 12.42 Å ³ ( ρ = 0.999 g cm −3) for classical and PI molecular-dynamics simulations, respectively. Path-integral simulations were performed using 32 beads per atoms. The water surface was simulated with a slab containing 108-molecules, in a cell of dimensions 12.74 × 12.74 × 57.35 Å ³. After equilibration, we carried out simulations for 0.5-1 ns for the bulk and the slab. In particular, we insured that the electric field across the water slab vanished within 0.0005 eV/Å. We tested size effects on the position of energy levels by performing MD simulations for 256-molecule bulk water and 216- and 384-molecule slabs. The details of the simulations used for finite size-scaling are discussed in the Supplementary Methods. We note that the choice of the system size in our MD simulations was limited by the cost of MBPT calculations on MD trajectories.

Electronic structure calculations

We carried out calculations of the electronic properties of water and water slabs at the density-functional theory (DFT) and G ₀W ₀ level. DFT calculations were carried out using the Quantum Espresso code61, with a plane-wave energy cutoff of 85 Ry, and HSCV pseudopotentials62,63. G ₀W ₀ calculations were performed using the the West code64,65, starting from wavefunctions obtained with the PBE66, PBE067, sc-hybrid39, and RSH40 functionals. Quasiparticle energies were determined using the same exchange-correlation potential as the one employed in the self-consistent-field DFT calculations. All of our G ₀W ₀ calculations were carried out with 1600 eigenpotentials and extrapolated using the procedure described in the Supplementary Table 5. DFT calculations were done for 128 equally-spaced snapshots along each of the bulk trajectories. G ₀W ₀ calculations for bulk water, starting from hybrid functional wavefunctions, were performed for 4 snapshots. The G ₀W ₀ results for the slabs were obtained using computed DFT values for the slab and G ₀W ₀ corrections for the bulk. Our analysis of the accuracy of this approach and error estimates are given in the Supplementary Methods.

To determine absolute energy levels from DFT and G ₀W ₀ calculations performed in plane-wave basis sets, we aligned the plane-average electrostatic potentials in the bulk water and in the slab, and referenced them to the vacuum energy level, following the procedure proposed by van de Walle and Martin68, as explained in the Supplementary Methods.

Data availability

All data associated with our article are available at https://datahub.uchicago.edu/.

Additional information

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. 24.
Fang, C. et al. The accurate calculation of the band gap of liquid water by means of GW corrections applied to plane-wave density functional theory molecular dynamics simulations. Phys. Chem. Chem. Phys. 17, 365 (2015).

2. 35.
Gaiduk, A. P., Zhang, C., Gygi, F. & Galli, G. Structural and electronic properties of aqueous NaCl solutions from ab initio molecular dynamics simulations with hybrid density functionals. Chem. Phys. Lett. 604, 89 (2014).

3. 36.
Kurahashi, N. et al. Photoelectron spectroscopy of aqueous solutions: streaming potentials of NaX (X = Cl, Br, and I) solutions and electron binding energies of liquid water and X −. J. Chem. Phys. 140, 174506 (2014).

4. 41.
Gaiduk, A. P., Paesani, F. & Galli, G. Electron affinities of water clusters from density-functional and many-body-perturbation theory, arXiv e-Prints, http://arxiv.org/abs/1711.10948 (2017).

5. 44.
Wan, Q., Spanu, L., Galli, G. A. & Gygi, F. Raman spectra of liquid water from ab initio molecular dynamics: vibrational signatures of charge fluctuations in the hydrogen bond network. J. Chem. Theory Comput. 9, 4124 (2013).

6. 55.
Coe, J. V. et al. Photoelectron spectroscopy of hydrated electron cluster anions, (H ₂O)- n=2-69. J. Chem. Phys. 92, 3980 (1990).

7. 57.
Zhang, L., Zhu, D., Nathanson, G. M. & Hamers, R. J. Selective Photoelectrochemical Reduction of Aqueous CO ₂ to CO by Solvated Electrons. Angew. Chem. Int. Ed. 53, 9746 (2014).

8. 58.
Elles, C. G., Jailaubekov, A. E., Crowell, R. A. & Bradforth, S. E. Excitation-energy dependence of the mechanism for two-photon ionization of liquid H ₂O and D ₂O from 8.3 to 12.4 eV. J. Chem. Phys. 125, 044515 (2006).

9. 65.
West Code (version 1.1.0): http://west-code.org/ (retrieved Tue Dec 19 18:34:47 2017).

Download references

Acknowledgements

The authors gratefully acknowledge helpful discussions with Stephen Bradforth, Nicholas Brawand, Federico Giberti, François Gygi, Ryan McAvoy, and Robert Seidel. A.P.G., M.G., and G.G. were supported by MICCoM as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), Materials Sciences and Engineering Division (5J-30161-0010A). A.P.G. was also supported by the postdoctoral fellowship from the Natural Sciences and Engineering Research Council of Canada. T.A.P. was supported by the Lawrence Fellowship. F.P. was supported by the National Science Foundation through grant CHE-1453204 and used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation through grant ACI-1053575. Part of this work was performed under the auspices of the U.S. DOE at Lawrence Livermore National Laboratory under contract DE-AC52-07A27344. An award of computer time was provided by the INCITE program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DEAC02-06CH11357.

Author information

Affiliations

1. Institute for Molecular Engineering, The University of Chicago, Chicago, IL, 60637, USA

Alex P. Gaiduk

, Marco Govoni

& Giulia Galli

2. Lawrence Livermore National Laboratory, Livermore, CA, 94551, USA

3. Materials Science Division, Argonne National Laboratory, Argonne, IL, 60439, USA

Marco Govoni

& Giulia Galli

4. Department of Chemistry and Biochemistry, Materials Science and Engineering, San Diego Supercomputer Center, University of California, San Diego, 92093, USA

Contributions

A.P.G., T.A.P., F.P., and G.G. designed the research. T.A.P. and F.P. performed MB-pol simulations. A.P.G. performed most of the electronic-structure calculations and data analysis. M.G. implemented the G ₀W ₀ method and DDH functionals in the West code. A.P.G. and G.G. wrote the manuscript with contributions from F.P. and T.A.P. All authors contributed to the discussion of the results.

Competing interests

The authors declare no competing financial interests.

Corresponding authors

Correspondence to Francesco Paesani or Giulia Galli.

Electronic supplementary material

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

# data[].attributes.attachment_link.description
The electron affinity of liquid water is a fundamental property which has not yet been accurately measured. Here, the authors predict this property by coupling path-integral molecular dynamics with ab initio potentials and electronic structure calculations, revisiting several estimates used in the literature.