Explicit treatment of hydrogen bonds in the Universal Force Field: validation and application for Metal-Organic Frameworks, hydrates and host-guest complexes

TL;DR

This paper introduces an explicit hydrogen bond treatment within the Universal Force Field, improving accuracy in modeling hydrogen-bonded systems like metal-organic frameworks and hydrates, and enabling faster, more reliable computational analyses.

Contribution

The authors develop and validate a new explicit hydrogen bond approach in UFF, enhancing structural prediction accuracy for complex hydrogen-bonded systems without arbitrary charge assignment.

Findings

Reduces maximum error in cell parameters from 66% to 14%.

Decreases mean unsigned error from 14% to 4%.

Enables rapid and accurate modeling of hydrogen-bonded frameworks.

Abstract

A straightforward means to include explicit hydrogen bonds within the Universal Force Field is presented. Instead of treating hydrogen bonds as non-bonded interaction subjected to electrostatic and Lennard-Jones potentials, we introduce an explicit bond with negligible bond order, thus maintaining the structural integrity of the H-bonded complexes and avoiding the necessity to assign arbitrary charges to the system. The explicit hydrogen bond changes the coordination number of the acceptor site and the approach is thus most suitable for systems with under-coordinated atoms, such as many metal-organic frameworks, however, it also shows excellent performance for other systems involving a hydrogen-bonded framework. In particular, it is an excellent means for creating starting structures for molecular dynamics and for investigations employing more sophisticated methods. The approach is validated for the hydrogen bonded complexes in the S22 dataset and then employed for a set of metal-organic frameworks from the Computation-Ready Experimental (CoRE) database and several hydrogen bonded crystals including water ice and clathrates. We show that direct inclusion of hydrogen bonds reduces the maximum error in predicted cell parameters from 66% to only 14% and the mean unsigned error is similarly reduced from 14% to only 4%. We posit that with the inclusion of hydrogen bonding, the solvent-mediated breathing of frameworks such as MIL-53 is now accessible to rapid UFF calculations, which will further the aim of rapid computational scanning of metal-organic frameworks while providing better starting points for electronic structure calculations.

Tables

Table 1. Table 1:  Hydrogen bond distances calculated using UFF for hydrogen bonded complexes in the S22 database Jurecka  et al. ( 2006 ) . Percent errors are calculated as ( X UFF subscript 𝑋 UFF X_{\textrm{UFF}} - X ref subscript 𝑋 ref X_{\textrm{ref}} )/ X ref × 100 subscript 𝑋 ref 100 X_{\textrm{ref}}\times 100 , where X UFF subscript 𝑋 UFF X_{\textrm{UFF}} denotes UFF-predicted and X ref subscript 𝑋 ref X_{\textrm{ref}} denotes the original ab initio value.

Complex (symmetry)	\pbox4cmHydrogen bond
Atom types	\pbox3cmUFF with
HF/3-21G
charges	\pbox3cmUFF with
explicit
H-bonds	Reference	% error
\ce(H2O)2 ($C_s$)	O_3-H$\cdots$O_3	2.666	1.920	1.952	-1.7
ammonia dimer ($C_{2h}$)	N_3-H$\cdots$N_3	2.854	-	2.504	-
Formic acid dimer ($C_{2h}$)	O_R-H$\cdots$O_1	2.488	1.886	1.670	12.9
Formamide dimer ($C_{2h}$)	N_R-H$\cdots$O_1	2.520	1.884	1.840	2.4
Uracil dimer ($C_{2h}$)	N_R-H$\cdots$O_2	2.481	1.882	1.774	6.1
2-pyridoxine - 2-aminopyridine ($C_1$)	N_R-H$\cdots$N_R	2.609	2.030	1.860	9.1
	N_R-H$\cdots$O_1	2.585	1.858	1.874	-0.8
Adenine - thymine WC ($C_1$)	N_R-H$\cdots$N_R	2.556	2.030	1.819	11.6
	N_R-H$\cdots$O_1	2.477	1.857	1.929	-3.7

Table 2. Table 2:  Error on bond angles calculated using UFF for hydrogen bonded complexes in the S22 database Jurecka  et al. ( 2006 ) . Percent errors are calculated as 100 / N ang × ∑ i = 1 N ang ( X UFF − X ref ) / X ref 100 subscript 𝑁 ang superscript subscript 𝑖 1 subscript 𝑁 ang subscript 𝑋 UFF subscript 𝑋 ref subscript 𝑋 ref 100/N_{\textrm{ang}}\times\sum\limits_{i=1}^{N_{\textrm{ang}}}{(X_{\textrm{UFF}}-X_{\textrm{ref}})/X_{\textrm{ref}}} and rmse as ( 1 / N ang × ∑ i = 1 N ang ( X UFF − X ref ) 2 ) 1 2 superscript 1 subscript 𝑁 ang superscript subscript 𝑖 1 subscript 𝑁 ang superscript subscript 𝑋 UFF subscript 𝑋 ref 2 1 2 (1/N_{\textrm{ang}}\times\sum\limits_{i=1}^{N_{\textrm{ang}}}{(X_{\textrm{UFF}}-X_{\textrm{ref}})^{2}})^{\frac{1}{2}} , where X UFF subscript 𝑋 UFF X_{\textrm{UFF}} denotes UFF-predicted and X ref subscript 𝑋 ref X_{\textrm{ref}} denotes the original ab initio value

	UFF, LJ only		UFF with explicit H-bonds
Complex (symmetry)	RMSE	Mean error (%)	RMSE	Mean error (%)
\ce(H2O)2 ($C_s$)	2.97341	2.6813	3.60928	2.85337
ammonia dimer ($C_{2h}$)	6.89258	7.00757	22.4794	21.6555
Formic acid dimer ($C_{2h}$)	6.97037	22.1015	7.08619	81.4685
Formamide dimer ($C_{2h}$)	9.29674	39.9108	2.85878	14.4514
Uracil dimer ($C_{2h}$)	3.60883	9.89631	3.02618	9.90769
2-pyridoxine - 2-aminopyridine ($C_1$)	3.57988	16.3517	2.85318	10.1884
Adenine - thymine WC ($C_1$)	3.7847	15.6545	4.11144	11.9441

Table 3. Table 3: Comparison of UFF calculated and experimental cell parameters of selected MOFs. The first 10 structures possess primarily intermolecular hydrogen bonds, the second 10 possess more intramolecular hydrogen bonds. Percent errors are calculated as ( X UFF subscript 𝑋 UFF X_{\textrm{UFF}} - X exp subscript 𝑋 exp X_{\textrm{exp}} )/ X exp × 100 subscript 𝑋 exp 100 X_{\textrm{exp}}\times 100 , where X UFF subscript 𝑋 UFF X_{\textrm{UFF}} denotes UFF-predicted and X exp subscript 𝑋 exp X_{\textrm{exp}} denotes the original value.

CSD Refcode		Experimental	\pbox3cmUFF without
H-bonds	\pbox3cmUFF with
H-bonds	\pbox3cm % error without
H-bonds	\pbox3cm% error with
H-bonds
MIL-53(Al) NPLoiseau et al. (2004)	a =	19.504	18.593	19.423	-4.7	-0.4
	b =	15.201	21.839	15.427	43.7	1.5
	c =	6.569	6.366	6.470	-3.1	-1.5
CDLGLU01Tomlin et al. (1997)	a =	11.575	10.099	11.217	-12.8	-3.1
	b =	10.764	15.304	9.920	42.2	-7.8
	c =	7.256	7.435	7.775	2.5	7.2
CUGLTMGramaccioli and Marsh (1966)	a =	11.084	13.404	10.989	20.9	-0.9
	b =	10.350	11.262	9.412	8.8	-9.1
	c =	7.238	6.330	7.246	-12.6	0.1
FUFREELi et al. (2009)	a =	7.785	8.312	7.808	6.8	0.3
	b =	10.238	12.172	9.780	18.9	-4.5
	c =	15.851	15.524	16.074	-2.1	1.4
HURNOXZhang et al. (2003)	a =	7.149	6.891	6.158	-3.6	-13.9
	b =	10.468	14.307	9.573	36.7	-8.5
	c =	11.295	10.156	11.518	-10.1	2.0
ICOWONBromant et al. (2006)	a =	8.754	8.240	8.232	-5.9	-6.0
	b =	10.003	10.693	9.775	6.9	-2.3
	c =	11.790	11.446	11.479	-2.9	-2.6
JUCXEKGrigorev et al. (1991)	a =	11.382	13.971	12.306	22.7	8.1
	b =	11.382	14.020	12.509	23.2	9.9
	c =	11.734	9.592	10.169	-18.3	-13.3
MAZTIRGuo and Guo (2006)	a =	19.014	19.839	19.215	4.3	1.1
	b =	6.843	9.715	6.782	42.0	-0.9
	c =	9.420	9.217	9.753	-2.2	3.5
MEHPAQMizutani et al. (1998)	a =	7.221	6.350	7.205	-12.1	-0.2
	b =	10.295	11.218	9.948	9.0	-3.4
	c =	11.041	13.425	10.966	21.6	-0.7
MUTVUTHuang et al. (2010)	a =	9.335	9.255	9.399	-0.9	0.7
	b =	22.203	25.664	23.168	15.6	4.3
	c =	27.514	26.756	27.710	-2.8	0.7
SALLATWu and Shi (2004)	a =	15.686	12.126	15.023	-22.7	-4.2
	b =	8.165	13.522	8.865	65.6	8.6
	c =	13.119	12.225	12.767	-6.8	-2.7
DUQSEORodriguez-Dieguez et al. (2010)	a =	8.884	7.930	7.606	-10.7	-14.4
	b =	13.093	11.963	12.266	-8.6	-6.3
	c =	13.135	16.660	14.021	26.8	6.7
FAPTUNWang et al. (2012a)	a =	14.304	11.704	14.757	-18.2	3.2
	b =	16.970	17.091	16.907	0.7	-0.4
	c =	11.098	11.342	11.373	2.2	2.5
KOJCUILi et al. (2008)	a =	18.326	21.776	18.356	18.8	0.2
	b =	25.300	21.273	24.369	-15.9	-3.7
	c =	7.524	6.804	6.900	-9.6	-8.3
RATVEPZheng et al. (2012)	a =	7.764	9.142	8.044	17.8	3.6
	b =	10.177	11.938	10.240	17.3	0.6
	c =	15.949	15.840	16.222	-0.7	1.7
SARBOEWang et al. (2012b)	a =	14.389	14.925	14.477	3.7	0.6
	b =	15.475	14.511	15.006	-6.2	-3.0
	c =	8.424	10.610	8.788	26.0	4.3
SIVKAKKawasaki et al. (2008)	a =	7.583	9.243	8.340	21.9	10.0
	b =	31.086	27.033	29.315	-13.0	-5.7
	c =	13.799	14.199	13.523	2.9	-2.0
VEFLUPXu and Xu (2012)	a =	8.002	7.283	7.460	-9.0	-6.8
	b =	9.091	10.877	9.031	19.6	-0.7
	c =	12.935	13.070	13.071	1.0	1.0
YORZAHDas et al. (2009)	a =	16.493	16.283	16.334	-1.3	-1.0
	b =	6.820	7.905	6.597	15.9	-3.3
	c =	22.220	23.331	21.906	5.0	-1.4
YUVSUEAn et al. (2010)	a =	15.435	17.188	15.518	11.4	0.5
	b =	15.435	17.188	15.793	11.4	2.3
	c =	22.775	18.122	21.005	-20.4	-7.8
ZNGLUDGramaccioli (1966)	a =	11.190	13.840	12.126	23.7	8.4
	b =	10.463	11.273	10.112	7.7	-3.4
	c =	7.220	6.369	7.203	-11.8	-0.2
				Maximum Unsigned Error	65.6	14.4
				Mean Unsigned Error	13.9	4.0

Table 4. Table 4:  Comparison of UFF calculated and experimental cell parameters of selected hydrogen bonded crystals. Percent errors are calculated as ( X UFF subscript 𝑋 UFF X_{\textrm{UFF}} - X exp subscript 𝑋 exp X_{\textrm{exp}} )/ X exp × 100 subscript 𝑋 exp 100 X_{\textrm{exp}}\times 100 , where X UFF subscript 𝑋 UFF X_{\textrm{UFF}} denotes UFF-predicted and X exp subscript 𝑋 exp X_{\textrm{exp}} denotes the original value.

CSD Refcode		Experimental	\pbox3cmUFF without
H-bonds	\pbox3cmUFF with
H-bonds	\pbox3cm % error without
H-bonds	\pbox3cm% error with
H-bonds
Ic iceEngel et al. (2015)	a	6.358	8.805	6.839	38.5	7.6
	b	6.358	9.019	6.641	41.8	4.5
	c	6.358	9.019	6.641	41.8	4.5
Ih iceEngel et al. (2015)	a	4.506	5.483	4.756	21.7	5.5
	b	4.506	6.284	4.717	39.5	4.7
	c	7.346	5.205	7.776	-29.2	5.9
alpha methanolTorrie et al. (1989)	a	4.873	4.645	4.955	-4.7	1.7
Figure 2(c)	b	4.641	5.000	4.487	7.7	-3.3
	c	8.867	9.835	7.798	10.9	-12.1
alpha methanolTorrie et al. (1989)				4.876		0.1
Figure 2(b)				5.178		11.6
				8.592		-3.1
beta methanolTauer and Lipscomb (1952)	a	6.409	4.362	6.702	-31.9	4.6
	b	7.199	7.268	7.099	1.0	-1.4
	c	4.649	6.895	4.554	48.3	-2.0
ammoniaOlovsson and Templeton (1959)	a	5.138	5.313	4.919	3.4	-4.3
	b	5.138	5.313	4.919	3.4	-4.3
	c	5.138	5.313	4.919	3.4	-4.3
Methane Hydrate-ILenz and Ojamäe (2011)	a	11.620	12.181	12.492	4.8	7.5
	b	11.620	16.164	12.490	39.1	7.5
	c	11.620	10.436	12.431	-10.2	7.0
Methane Hydrate-IILenz and Ojamäe (2011)	a	11.890	12.262	12.745	3.1	7.2
	b	11.890	11.775	12.723	-1.0	7.0
	c	11.890	15.180	12.712	27.7	6.9
Methane Hydrate-HLenz and Ojamäe (2011)	a	11.910	10.826	12.789	-9.1	7.4
	b	11.910	14.074	12.813	18.2	7.6
	c	9.894	11.302	10.456	14.2	5.7
WUVZIWLindner et al. (2003)	a	7.400	7.393	7.355	-0.1	-0.6
	b	24.448	24.958	25.248	2.1	3.3
	c	14.265	15.399	14.707	8.0	3.1
PgC2 Cocrystal 5 Fowler et al. (2014)	a	10.059	9.642	9.693	-4.1	-3.6
	b	14.556	15.927	13.886	9.4	-4.6
	c	15.558	17.731	14.945	14.0	-3.9
PgC4 Cocrystal 6Fowler et al. (2014)	a	21.075	23.775	21.334	12.8	1.2
	b	12.796	13.363	12.452	4.4	-2.7
	c	18.229	18.353	17.438	0.7	-4.3
PgC4 Cocrystal 7Fowler et al. (2014)	a	21.075	23.490	20.088	11.5	-4.7
	b	12.796	13.506	12.380	5.6	-3.2
	c	18.229	18.251	17.118	0.1	-6.1
ABAZOSYeo and Harris (1999)	a	16.338	17.022	16.306	4.2	-0.2
	b	21.933	22.805	20.839	4.0	-5.0
	c	16.338	16.870	16.310	3.3	-0.2
ABEBUFGordon-Wylie et al. (2004)	a	10.722	10.989	10.851	2.5	1.2
	b	10.900	11.831	11.230	8.5	3.0
	c	27.635	25.745	27.648	-6.8	0.0
ABUCIJApel et al. (2001)	a	11.281	10.924	11.076	-3.2	-1.8
	b	17.888	18.370	17.973	2.7	0.5
	c	23.950	24.319	23.863	1.5	-0.4
SIYRAUTelfer et al. (2008)	a	13.456	13.861	13.257	3.0	-1.5
	b	14.394	13.727	14.647	-4.6	1.8
	c	15.935	18.019	16.945	13.1	6.3
WARWOBChang et al. (1993)	a	19.297	22.634	20.706	17.3	7.3
	b	4.616	5.599	4.589	21.3	-0.6
	c	8.705	6.954	7.972	-20.1	-8.4
				Maximum Unsigned Error	48.3	12.1
				Mean Unsigned Error	12.6	4.2

Equations

$$r_{ij}=r_i+r_j+r_{\text{BO}}-r_{\text{EN}}$$

$$r_{\text{BO}}=-\lambda (r_i+r_j)\ln (n)$$

Full text

The approach is validated for the hydrogen bonded complexes in the S22 dataset and then employed for a set of metal-organic frameworks from the Computation-Ready Experimental (CoRE) database and several hydrogen bonded crystals including water ice and clathrates. We show that direct inclusion of hydrogen bonds reduces the maximum error in predicted cell parameters from 66% to only 14% and the mean unsigned error is similarly reduced from 14% to only 4%. We posit that with the inclusion of hydrogen bonding, the solvent-mediated breathing of frameworks such as MIL-53 is now accessible to rapid UFF calculations, which will further the aim of rapid computational scanning of metal-organic frameworks while providing better starting points for electronic structure calculations.

I Introduction

Metal-Organic Frameworks (MOFs) are now a well-known class of crystalline, porous materials, where inorganic connectors are joined by organic linkers, forming a 3-dimensional network. Since the publication of MOF-5 in 1999Li et al. (1999), research into MOFs has grown almost exponentially and there are now several thousand synthesised MOFs.Chung et al. (2014) Several reviews have been written on various aspects of MOF chemistry,Odoh et al. (2015); Schneemann et al. (2014) including their interactions with water.Burtch et al. (2014); Canivet et al. (2014)

Water is an important solvent for MOFs, it serves as a synthetic medium,Choi et al. (2004); Yang et al. (2007); Stock and Biswas (2012) a structural component,DeCoste et al. (2013) and an adsorbate.Furukawa et al. (2014) Water plays a major role in mediating defects and dissolution of many MOFs.St. Petkov et al. (2012); Tan et al. (2015); Küsgens et al. (2009) Some MOFs, such as the UiO series, are highly water stableCavka et al. (2008), while others break down or transformSaha and Deng (2010); Greathouse and Allendorf (2006); Cheng et al. (2009); Bezverkhyy et al. (2014) when exposed to water. Proposed applications of MOFs include separation of organic contaminantsWang et al. (2014a) and arsenicVu et al. (2015) from aqueous mixtures and capture of water from airFurukawa et al. (2014).

The requirement to computationally scan or screen the nearly infinite number of possible MOFs led to the development of UFF4MOF.Addicoat et al. (2014); Coupry et al. (2016) UFF4MOF extends the Universal Force Field of RappéRappe et al. (1992) to include several atom types present in MOFs but not accounted for in the original list of atom types. UFF4MOF thus permits the rapid calculation of structures for a wide range of actual and hypothetical MOFs. The structures predicted using UFF4MOF are typically within 5% of experimental cell parameters, and it is therefore of interest whether similar accuracy can be maintained when the structural parameters are significantly determined by hydrogen bonding within the framework. MIL-53Loiseau et al. (2004) represents a notable case of water not only stabilizing, but determining the framework structure. The transition from (evacuated) large pore to narrow pore to hydrated and superhydrated large pore has been extensively studied by both experimental and computational methods.Salles et al. (2011); Cirera et al. (2012); Guillou et al. (2011); Ortiz et al. (2014) The empty, large pore (LP) structures of MIL-53(M) (M = Cr, Al, Fe) have all been replicated using UFF4MOF, and in this work we seek to replicate the hydrated narrow pore (NP) structures, whose cell parameters are dictated by the hydrogen bonds that occur between the guest water molecules and the framework.Salazar et al. (2015); Ferey (2012)

Treatment of hydrogen bonds is challenging for classical force fields. On one hand, the bonds are readily formed and broken in solution, so that an explicit, permanent bond would be an inappropriate description if investigating dynamical processes or, for example, the liquid state. Many force fields treat hydrogen bonds as non-bonded interactions: the bonding is achieved by attractive components arising from the London dispersion and the Coulomb interaction. In UFF, the former is defined by a Lennard-Jones potential, while the latter relies on the definition of charges. As the Coulomb interaction dominates, charge definition is crucial for the results. However, this approach works well for water and is particularly useful in QM/MM calculations, as was previously shown in DFTB/UFF calculations with electrostatic embedding.Heine et al. (2007); Lula et al. (2007).

Several force fields attempt to explicitly treat hydrogen bonds: MM2Allinger (1977), MM3Allinger et al. (1989) and AMBER force fields all employ (optionally in the case of AMBER) an angle-independent 10,12 Lennard-Jones potential to describe nonbonded interactions. MM3 was later updated to include directional hydrogen bondingLii and Allinger (1994). Paton and GoodmanPaton and Goodman (2009) reviewed these force fields along with OPLS*, OPLSAA and MMFF, which do not include explicit hydrogen bonding and concluded that for the S22 database, the OPLSAA and MMFF force fields performed best, yielding accurate geometries and interaction energies. The same two force fields also ranked highest against the JSCH2005 database, but some optimized geometries differed from the benchmark geometry.

Despite these advances, there generally remains a technical problem. For host-guest systems that are treated entirely using a force field, it is often not practical to define charges to treat the electrostatic interactions. Indeed, although a charge equilibration scheme was designed for use with UFF Rappe and Goddard (1991), the original parameterization was done without an electrostatic model: the implementation is thus open to interpretation. Some codes forego explicit charges completely Martin (2013); O’Boyle et al. (2011). UFF also explicitly includes electrostatic interactions for 1st and 2nd neighbours, so that most electrostatics is screened for a bonded system. This, of course, does not hold for non-bonded interactions, where significant electrostatic contributions arise due to the proximity of atoms formally carrying a charge. There is a further, merely technical aspect for the choice of explicit hydrogen bond that is crucial when studying static hydrogen-bonded frameworks including MOFs, ice, clathrates, or inclusion compounds: upon generation their starting structure, the non-bonded character in conjunction with substantial forces at initial geometry optimization, or, worse, kinetic energy gain upon molecular dynamics startup, makes it hard to converge to or to maintain the anticipated local minimum representing the desired framework topology.

II Definition of explicit hydrogen bonds within the Universal Force Field

To provide a fast and topology-preserving structure generation method, and to avoid the definition of charges in order to account for the Coulomb interaction in the hydrogen bonds, we take advantage of the implicit definition of the electronegativity correction $r_{\textrm{EN}}$, which is contributing to the UFF natural bond length $r_{ij}$:

[TABLE]

where $r_{i}$ and $r_{j}$ are the bond radii of atoms $i$ and $j$ respectively, $r_{\textrm{BO}}$ is a bond order correction to the bond distance. The bond order correction, in turn, is defined as:

[TABLE]

where $\lambda$ is a proportionality constant derived using propane, propene and propyne with the C_3, C_2 and C_1 radii and is equal to 0.1332. There is precedent for assigning unique bond orders for special types of bonds, with the amide bond order set to 1.41 in order to reproduce the C-N bond distance in $N$-methylformamide.Rappe et al. (1992) Assigning a similar bond order for hydrogen bonds requires selecting an appropriate reference. At first glance, this is less than straightforward, given the variety of hydrogen bonds that may be encountered in framework structures. However, considering the prime importance of water to the synthesis, structure and behaviour of MOFs, a reasonable reference is the the prototypical $C_{s}$ global minimum of the water dimer. We further require that our treatment of hydrogen bonds can reproduce the hydrogen bond mediated breathing that occurs in hydrated MIL-53.

Treating hydrogen bonds as explicit bonds without bond order does have negative consequences: The dynamical breaking and formation of hydrogen bonded networks, essential for the description of the liquid state, becomes impossible (though it may work if the Reactive Force Field, ReaxFFvan Duin et al. (2001), was used instead of UFF). We note that defining an explicit hydrogen bond increases the coordination number of the central atom by one, however, this typically does not pose a problem for otherwise under-coordinated sites, and fully coordinated atoms typically are not hosting extra solvent. Secondly, because the defined bonding network is changed in the hydrogen-bonded system, it is not possible to calculate binding energies or relative stabilities using this approach. However, for the problem of rapid structural pre-optimization, neither of these disadvantages apply and furthermore, they are offset by the ability to rapidly produce high quality structures, without requiring the relatively expensive calculation of electrostatic terms.

As we will show below, for “frozen” configurations of hydrogen bonded dimers our approach performs significantly better than the traditional UFF treatment using nonbonded interactions with Coulomb interactions defined via HF/3-21G charges, and should provide a much better starting point for subsequent electronic structure calculations.

The water dimer has been extensively studied over decades, by both experimentalBouteiller et al. (2011); Serov et al. (2014) and computationalSong et al. (2014); Wang et al. (2014b); Ronca et al. (2014) means. Early calculations yielded a H$\cdots$O distance of 1.72 ÅMorokuma and Pedersen (1968), a subsequent study including some of the same authors yielded a value of 1.80 ÅMorokuma and Winick (1970). By the 1990’s, calculations using MP2 and Coupled Electron Pair Theory (CEPA-1)van Duijneveldt-van de Rijdt and van Duijneveldt (1992) yielded geometries and energies in good agreement with experimentOdutola and Dyke (1980); Dyke et al. (1977) and the “gold standard” CCSD(T)/QZ geometry found in the S22 databaseJurecka et al. (2006) also agrees very well, with $R_{H\cdots O}=1.952\AA$ (corresponding to $R_{OO}=2.91\AA$). If we therefore consider the UFF bond length and bond order correction in Equations 1 and 2, we note that $\lambda$, $r_{i}$ and $r_{j}$ and $r_{\textrm{EN}}$ are all fixed within the UFF framework. Therefore, for water, where the oxygen atom has the O_3 atom type and hydrogen has the H_ atom type, $r_{i}+r_{j}=0.354+0.658=1.012\AA$, the electronegativity correction for a H-O_3 bond is 0.0021 and thus the required bond order correction is equal to 0.9380 Å. A bond order of 0.001 yields correction of 0.9312 Å and noting the diminishing returns of further optimization of such a “bond order”, combined with the large range of acceptable bond lengths for a hydrogen bond, we propose to this value as a reasonable approximation to a bond order for describing hydrogen bonds in UFF calculations on metal-organic frameworks. Employing this bond order yields a bond length of 2.025 Å for a H$\cdots$N_R bond and 1.899 Å where an O_2 atom is the proton acceptor, both of which are reasonable lengths for a hydrogen bond.

III Results and Discussion

To check the basic sensibility of using a bond order to correct hydrogen bond distances, we undertook geometry optimizations of the hydrogen bonded complexes in the S22 database. The $C_{2h}$ ammonia dimer was excluded from analysis as the angle terms resulting from the addition of the two hydrogen bonds considerably change the geometry of the dimer. The hydrogen bond distances of the other six complexes are shown in Table 1 and this is the only metric we employ for these non-bonded clusters.Witte et al. (2015) Calculations were undertaken in the General Utility Lattice Program (GULP)Gale and Rohl (2003); Gale (2005) except those employing atomic charges, where deMonNanodeM (2009) was employed. Hartree-Fock calculations were undertaken in Gaussian09Frisch et al. (2009).

For these simple complexes, the results show reasonable agreement with the ab initio reference and are far superior to UFF + Lennard-Jones + electrostatics. The two N_R-H$\cdots$N_R bonds are overestimated by approximately 10%, as is the unusually short hydrogen bond in the formic acid dimer. Other bonds are within 6% of their reference values. Undertaking the calculations without specification of the hydrogen bond results in bond distances increasing by approximately 1Å  which is clearly poor.

We note in Table 1, the only failure of the explicit hydrogen-bond approach is the $C_{2h}$ ammonia dimer, which upon inclusion of two explicit hydrogen bonds, optimizes to a singly hydrogen-bonded complex of $C_{s}$ symmetry. In this case the two hydrogen bonds in the reference structure form H-N$\cdots$H angles of 58∘, a significant deviation from the 106.7∘ angle of the N_3 parameter, and the optimizer prefers to allow a single hydrogen bond with a close-to-ideal angle (i.e. essentially tetrahedral geometry around the acceptor nitrogen atom), rather than two hydrogen bonds with large errors on the angle term. This is a general limitation of this approach, in that by making the hydrogen bond explicit, the coordination number of both the hydrogen atom and the accepting atom are increased by one and the hydrogen bonded atom consequently figures in the angle terms around the acceptor atom. In most cases, this effect is either desired or benign, such as in the case of hydrogen bonding to an under-coordinated metal atom in a paddlewheel.

III.1 Metal-Organic Frameworks

Having thus established the validity of the approach, a set of framework materials where hydrogen bonding is important for maintaining structural integrity was selected from the Computation Ready Experimental (CoRE) database.Chung et al. (2014) The original crystal structures were re-sourced from the Cambridge Structural DatabaseAllen (2002) in order to recover the solvent molecules. To this test set we add the particular case of MIL-53(Al) NPLoiseau et al. (2004). After assigning atom types to each structure, we detect hydrogen bonds by looking for hydrogen atoms. For every H atom that is directly connected to (O, S, N) we assign a larger covalent radius (of 2Å), then with this, we re-build the connectivity list and we add to the original connectivity list every new bond between these super-big hydrogen atoms and (O, N, S, F, Cl, Br, I) provided that the bond angle is between 140 and 220 degrees. For all structures we calculate the structure both with and without specification of hydrogen bonds and note the resultant cell parameters, which are listed in Table LABEL:tab:CoRE_hbonds.

In framework systems, the explicit inclusion of hydrogen bonds should yield even better agreement with experimental reference structure, as terms corresponding to regular covalent bonds far outnumber the contribution due to hydrogen bonds, and this is indeed the case. Without specifying hydrogen bonds, each structure has at least one cell parameter expand by 15% or greater, and the mean unsigned error on all cell parameters is 13.9%. Once hydrogen bonds are specified, the maximum unsigned error in cell parameters is only 14.4% and the mean unsigned error is only 4.0%.

The narrow pore structure of MIL-53(Al) is of particular note. Each pore contains two water molecules which are hydrogen-bonded to the framework oxygen atoms and the hydrogen of the framework hydroxy group. Optimising the structure with 16 hydrogen bonds specified (two hydrogen bonds per water molecule, see Figure 1), the cell parameters of MIL-53(Al) are all predicted within 1.5%, including the $b$ dimension, which otherwise expands by over 43% to resemble the large pore structure.

III.2 Hydrogen-bonded crystal structures

The prediction of the structures of hydrogen-bonded MOFs is perhaps the primary use-case of this approach, however, the structure of any static hydrogen-bonded system is amenable to calculation this way. To illustrate the approach’ broader utility, a variety of structures where hydrogen bonding is crucial to defining the structure were recalculated. Hydrogen bonds were inserted following the literature description of each structure. Table 4 shows the results for these structures.

III.2.1 Methanol

The case of $\alpha$-methanolTorrie et al. (1989) deserves special consideration as it illustrates how the choice of hydrogen bonds affects the optimized structure. There are two distinct types of hydrogen bond in $\alpha$-methanol; a bond from the hydroxyl hydrogen to the neighbouring oxygen with a OH$\cdots$O distance of 1.75Å and a CH$\cdots$O distance of 2.51Å, both illustrated in Figure 2. Choosing only the strong OH$\cdots$O results in a correct estimate of the $a$ parameter, but exchanges the errors on the $b$ and $c$ parameters. Including the CH$\cdots$O in addition results in an underestimation of both the $b$ and $c$ parameters and a moderate overestimation of $a$. While none of the three calculated structures is a clear winner, arguably including only the OH$\cdots$O results in the best quality structure. This choice is then consistent with the structure of $\beta$-methanolTauer and Lipscomb (1952), the high temperature phase, which possesses only OH$\cdots$O hydrogen bonds. These bonds, arranged in sheets, are crucial to the accurate description of the structure, neglecting them results in a strong distortion of the crystal, with the $a$ parameter being underestimated by 32% and the $c$ parameter being overestimated by 48%.

III.2.2 Water ice and hydrate structures

The classic hydrogen-bonded structure is that of water. The approach of fixing hydrogen bonds as actual bonds, makes the description of fluxional systems, such as liquid water, impossible. However, the structure of non-fluxional ice is readily amenable to calculation, the structures of both cubic and hexagonal water ice were calculated using all hydrogen bonds, and making each oxygen atom formally tetrahedral with two single bonds and two hydrogen bonds. Employing this approach results in a nearly uniform overestimation of lattice parameters of 4.5-7.5% for both structures, a clear improvement over neglecting the hydrogen bonds which results in parameters being over/underestimated by 20-40%.

The success of the explicit approach to hydrogen bonds in water ice has further significance in describing clathrate structures. Methane hydrates, in particular are the subject of renewed research as, found on the ocean floor and in polar regions they are estimated to contain up to 12% of all the organic carbon on EarthGbaruko et al. (2007), making them an important energy resource. In addition, hydrates are often formed within gas pipelines, where they are unwanted and cause significant damageSloan (2003). The structures of three methane hydrates, MH-I, MH-II and MH-H (hexagonal)Lenz and Ojamäe (2011) were calculated, yielding results broadly similar to those of water ice, whereby specification of hydrogen bonds gives a structure with cell parameters uniformly overestimated by approximately 7% and neglecting those bonds results in very poor structures with parameters under- and overestimated by up to 40%. More modest improvement is seen for the structure of $n$-butanol hexahydrate (CCSD Refcode WUVZIW).

III.2.3 Host-guest inclusion complexes and cocrystals

To further illustrate the diversity of hydrogen-bonded systems to which this simple approach may be applied, structures described as being hydrogen-bonded were sourced from the CCSDAllen (2002). Structures included two urea-based inclusion compounds (ABAZOSYeo and Harris (1999) and WARWOBChang et al. (1993)), a macrocyclic inclusion complex, ABUCIJApel et al. (2001), a porous diamide matrix, ABEBUFGordon-Wylie et al. (2004), a Cu coordination compound forming a hydrogen-bonded helicate, SIYRAUTelfer et al. (2008) and three cocrystals of pyrogallol[4]arenes and the ionic liquid 1-ethyl-3-methylimidazolium ethylsulfateFowler et al. (2014). In these cases, neglecting hydrogen bonds typically leads to at least one cell parameter in error by greater than 10%. Overall, employing explicit hydrogen bonds reduces the maximum unsigned error for these complexes from 48.3% to 12.1% and the mean unsigned error from 12.6% to 4.2%.

IV Potential Applications

When used to predict solvent effects, the accuracy of the structures produced by our approach is highly dependant on the guessed positions of the solvent molecules. For the structures where case by case reasoning is impossible, using a genetic algorithm (e.g: as implemented in the Atomic Simulation EnvironmentBahn and Jacobsen (2002); Vilhelmsen and Hammer (2014)) to maximize the number of hydrogen bonds in the structure of interest can yield reasonable starting points, enabling the use of explicit bonding in automated tasks. An active area of research like water confinement in carbon nanotubes of varying diameters Pascal et al. (2011), where the hydrogen bonding is the single most important factor for structure determination, would likely benefit from cheap, systematic structure elucidation. Appropriate modifications to the ASE code are underway to allow for solvent rigid motions during optimizations, and a proof of concept script is available on github (https://github.com/DCoupry/GenAlgHbond). However, the main application of explicit hydrogen bonding is case by case reasoning for a reasonable pre-optimization, followed by electronic structure methods. More complex properties of the structures, like frequencies (for which UFF was not parametrized) or the effects of H-bond anisotropy, fall squarely outside of the scope of this paper.

V Conclusions

We report a simple approach to explicitly treat hydrogen bonds within the Universal Force Field. The approach does not require any definition of atomic charges, which is a great technical (codes such as GULP and ADF do not include long-range and non-bonded electrostatic interactions within UFF) and practical (the definition of atomic charges is not well-defined) advantage. Moreover, it avoids the computationally costly electrostatic interaction term and is equally applicable for periodic and non-periodic systems.

The specification of hydrogen bonds with negligible (0.001) bond order increases the coordination number of the acceptor site by one, which is compatible with the definition of the main acceptor sites O and N, as well as with under-coordinated metal sites as present in metal-organic frameworks.

The approach was validated for the hydrogen bonded complexes in the S22 database. We show that it is very effective for both framework-framework hydrogen bonds and importantly, for framework-adsorbate bonds. Using this approach, both the large pore and narrow pore structures of MIL-53 can be calculated accurately with UFF. We show the generality of the approach by also applying it to hydrogen bonded crystal structures and host-guest inclusion complexes, including the environmentally and commercially important clathrates,

Finally, we note that this approach does not require any implementation, and thus works in any software that includes a UFF implementation. We propose this method to be most useful for pre-optimization and screening of static hydrogen-bonded systems.

VI Acknowledgements

This work was partially supported by the European Commission (ERC StG C3ENV (GA 256962), MSC-EID PROPAGATE (316897), DFG FOR 2433 FlexMOF (HE 3543/31-1) and the VolkswagenStiftung.