Towards improved FAIRness of the ThermoML Archive.

The ThermoML Archive is a subset of Thermodynamics Research Center (TRC) data holdings corresponding to cooperation between NIST TRC and five journals: Journal of Chemical Engineering and Data (ISSN: 1520-5134), The Journal of Chemical Thermodynamics (ISSN: 1096-3626), Fluid Phase Equilibria (ISSN: 0378-3812), Thermochimica Acta (ISSN: 0040-6031), and International Journal of Thermophysics (ISSN: 1572-9567). Data from initial cooperation (around 2003) through the 2019 calendar year are included. The archive has undergone a major update with the goal of improving the FAIRness and user experience of the data provided by the service. The web application provides comprehensive property browsing and searching capabilities; searching relies on a RESTful API provided by the Cordra software for managing digital objects. JSON files with a schema derived from ThermoML are provided as an additional serialization to lower the barrier to programmatic consumption of the information, for stakeholders who may have a preference of JSON over XML. The ThermoML and JSON files for all available entries can be downloaded from data.nist.gov (https://data.nist.gov/od/id/mds2-2422).

Survey of Data and Models for Refrigerant Mixtures Containing Halogenated Olefins.

We survey existing data for refrigerant blends containing halogenated olefins (hydrofluoroolefins (HFO), hydrochlorofluoroolefins (HCFO) and hydrochloroolefins (HCO)) in the open literature. The data are primarily taken from the NIST SOURCE database and are presented in graphical form to demonstrate the relative coverage of the data. The primary conclusion is that blends containing halogenated olefins remain only sparsely measured in experiments, and some classes of data (e.g., speed-of-sound data) are particularly sparse for blends containing halogenated olefins. The second part of this study compares the thermodynamic models in NIST REFPROP against the experimental datasets and identifies systems (of which there are many) where refitting of the thermodynamic model is required.

Environmental Mercury Chemistry - In Silico.

Mercury (Hg) is a global environmental contaminant. Major anthropogenic sources of Hg emission include gold mining and the burning of fossil fuels. Once deposited in aquatic environments, Hg can undergo redox reactions, form complexes with ligands, and adsorb onto particles. It can also be methylated by microorganisms. Mercury, especially its methylated form methylmercury, can be taken up by organisms, where it bioaccumulates and biomagnifies in the food chain, leading to detrimental effects on ecosystem and human health. In support of the recently enforced Minamata Convention on Mercury, a legally binding international convention aimed at reducing the anthropogenic emission of-and human exposure to-Hg, its global biogeochemical cycle must be understood. Thus, a detailed understanding of the molecular-level interactions of Hg is crucial. The ongoing rapid development of hardware and methods has brought computational chemistry to a point that it can usefully inform environmental science. This is particularly true for Hg, which is difficult to handle experimentally due to its ultratrace concentrations in the environment and its toxicity. The current account provides a synopsis of the application of computational chemistry to filling several major knowledge gaps in environmental Hg chemistry that have not been adequately addressed experimentally. Environmental Hg chemistry requires defining the factors that determine the relative affinities of different ligands for Hg species, as they are critical for understanding its speciation, transformation and bioaccumulation in the environment. Formation constants and the nature of bonding have been determined computationally for environmentally relevant Hg(II) complexes such as chlorides, hydroxides, sulfides and selenides, in various physical phases. Quantum chemistry has been used to determine the driving forces behind the speciation of Hg with hydrochalcogenide and halide ligands. Of particular importance is the detailed characterization of solvation effects. Indeed, the aqueous phase reverses trends in affinities found computationally in the gas phase. Computation has also been used to investigate complexes of methylmercury with (seleno)amino acids, providing a molecular-level understanding of the toxicological antagonism between Hg and selenium (Se). Furthermore, evidence is emerging that ice surfaces play an important role in Hg transport and transformation in polar and alpine regions. Therefore, the diffusion of Hg and its ions through an idealized ice surface has been characterized. Microorganisms are major players in environmental mercury cycling. Some methylate inorganic Hg species, whereas others demethylate methylmercury. Quantum chemistry has been used to investigate catalytic mechanisms of enzymatic Hg methylation and demethylation. The complex interplay between the myriad chemical reactions and transport properties both in and outside microbial cells determines net biogeochemical cycling. Prospects for scaling up molecular work to obtain a mechanistic understanding of Hg cycling with comprehensive multiscale biogeochemical modeling are also discussed.

Cys.sqlite: A Structured-Information Approach to the Comprehensive Analysis of Cysteine Disulfide Bonds in the Protein Databank.

Cysteine is a multifaceted amino acid that is central to the structure and function of many proteins. A disulfide bond formed between two cysteines restrains protein conformations through the strong covalent bond and torsions about the bond that prefer, energetically, ±90°. In this study, we transform over 30 000 Protein Databank files (PDBx/mmCIFs) into a single file, the SQLite database (Cys.sqlite). The database schema is designed to accommodate the structural information on both oxidized and reduced cysteines and to retain essential protein metadata to establish informational and biological provenance. Cys.sqlite contains over 95 000 peptide chains and 500 000 cysteines (700 000 structural conformers); there are over 265 000 cysteine disulfide bond conformations from structures solved with all available experimental methods. The structural information is analyzed with respect to sequence identity cutoff, the experimental method, and energetics of the disulfide. We find that as the experimental information becomes limiting and the influence of modeling becomes more pronounced, the observed average strain increases artificially. The database and analyses presented here can be used to improve the refinement of biological structures from experiments that are known to contain one or more disulfide bonds.

HackaMol: An Object-Oriented Modern Perl Library for Molecular Hacking on Multiple Scales.

HackaMol is an open source, object-oriented toolkit written in Modern Perl that organizes atoms within molecules and provides chemically intuitive attributes and methods. The library consists of two components: HackaMol, the core that contains classes for storing and manipulating molecular information, and HackaMol::X, the extensions that use the core. The core is well-tested, well-documented, and easy to install across computational platforms. The goal of the extensions is to provide a more flexible space for researchers to develop and share new methods. In this application note, we provide a description of the core classes and two extensions: HackaMol::X::Calculator, an abstract calculator that uses code references to generalize interfaces with external programs, and HackaMol::X::Vina, a structured class that provides an interface with the AutoDock Vina docking program.

X-ray structure of a Hg2+ complex of mercuric reductase (MerA) and quantum mechanical/molecular mechanical study of Hg2+ transfer between the C-terminal and buried catalytic site cysteine pairs.

Mercuric reductase, MerA, is a key enzyme in bacterial mercury resistance. This homodimeric enzyme captures and reduces toxic Hg2+ to Hg0, which is relatively unreactive and can exit the cell passively. Prior to reduction, the Hg2+ is transferred from a pair of cysteines (C558' and C559' using Tn501 numbering) at the C-terminus of one monomer to another pair of cysteines (C136 and C141) in the catalytic site of the other monomer. Here, we present the X-ray structure of the C-terminal Hg2+ complex of the C136A/C141A double mutant of the Tn501 MerA catalytic core and explore the molecular mechanism of this Hg transfer with quantum mechanical/molecular mechanical (QM/MM) calculations. The transfer is found to be nearly thermoneutral and to pass through a stable tricoordinated intermediate that is marginally less stable than the two end states. For the overall process, Hg2+ is always paired with at least two thiolates and thus is present at both the C-terminal and catalytic binding sites as a neutral complex. Prior to Hg2+ transfer, C141 is negatively charged. As Hg2+ is transferred into the catalytic site, a proton is transferred from C136 to C559' while C558' becomes negatively charged, resulting in the net transfer of a negative charge over a distance of ∼7.5 Å. Thus, the transport of this soft divalent cation is made energetically feasible by pairing a competition between multiple Cys thiols and/or thiolates for Hg2+ with a competition between the Hg2+ and protons for the thiolates.

Mercury methylation by HgcA: theory supports carbanion transfer to Hg(II).

Many proteins use corrinoid cofactors to facilitate methyl transfer reactions. Recently, a corrinoid protein, HgcA, has been shown to be required for the production of the neurotoxin methylmercury by anaerobic bacteria. A strictly conserved Cys residue in HgcA was predicted to be a lower-axial ligand to Co(III), which has never been observed in a corrinoid protein. Here, we use density functional theory to study homolytic and heterolytic Co-C bond dissociation and methyl transfer to Hg(II) substrates with model methylcobalamin complexes containing a lower-axial Cys or His ligand to cobalt, the latter of which is commonly found in other corrinoid proteins. We find that Cys thiolate coordination to Co facilitates both methyl radical and methyl carbanion transfer to Hg(II) substrates, but carbanion transfer is more favorable overall in the condensed phase. Thus, our findings are consistent with HgcA representing a new class of corrinoid protein capable of transferring methyl groups to electrophilic substrates.

Activation loop plasticity and active site coupling in the MAP kinase, ERK2.

Changes in the dynamics of the protein kinase, ERK2, have been shown to accompany its activation by dual phosphorylation. However, our knowledge about the conformational changes represented by these motions is incomplete. Previous NMR relaxation dispersion studies showed that active, dual-phosphorylated ERK2 undergoes global exchange between at least two energetically similar conformations. These findings, combined with measurements by hydrogen exchange mass spectrometry (HX-MS), suggested that the global conformational exchange involves motions of the activation loop (A-loop) that are coupled to regions surrounding the kinase active site. In order to better understand the contribution of dynamics to the activation of ERK2, we applied long conventional molecular dynamics (MD) simulations starting from crystal structures of active, phosphorylated (2P), and inactive, unphosphorylated (0P) ERK2. Individual trajectories were run for (5 to 25) µ s and totaled 727 µ s. The results showed that the A-loop is unexpectedly flexible in both 2P- and 0P-ERK2, and able to adopt multiple long-lived (>5 µ s) conformational states. Simulations starting from the X-ray structure of 2P-ERK2 (2ERK) revealed A-loop states corresponding to restrained dynamics within the N-lobe, including regions surrounding catalytic residues. One A-loop conformer forms lasting interactions with the C-terminal L16 segment and shows reduced RMSF and greater compaction in the active site. By contrast, simulations starting from the most common X-ray conformation of 0P-ERK2 (5UMO) reveal frequent excursions of A-loop residues away from a C-lobe docking site pocket and towards a new state that shows greater dynamics in the N-lobe and disorganization around the active site. Thus, the A-loop in ERK2 appears to switch between distinct conformational states that reflect allosteric coupling with the active site, likely occurring via the L16 segment. Analyses of crystal packing interactions across many structural datasets suggest that the A-loop observed in X-ray structures of ERK2 may be driven by lattice contacts and less representative of the solution structure. The novel conformational states identified by MD expand our understanding of ERK2 regulation, by linking the activated state of the kinase to reduced dynamics and greater compaction surrounding the catalytic site.

Activation Loop Plasticity and Active Site Coupling in the MAP Kinase, ERK2.

Previous studies of the protein kinase, ERK2, using NMR and hydrogen-exchange measurements have shown changes in dynamics accompanying its activation by phosphorylation. However, knowledge about the conformational motions involved is incomplete. Here, we examined ERK2 using long conventional molecular dynamics (MD) simulations starting from crystal structures of phosphorylated (2P) and unphosphorylated (0P) forms. Individual trajectories were run for (5 to 25) µs, totaling 727 µs. The results show unexpected flexibility of the A-loop, with multiple long-lived (>5 µs) conformational states in both 2P- and 0P-ERK2. Differential contact network and principal component analyses reveal coupling between the A-loop fold and active site dynamics, with evidence for conformational selection in the kinase core of 2P-ERK2 but not 0P-ERK2. Simulations of 2P-ERK2 show A-loop states corresponding to restrained dynamics within the N-lobe, including regions around catalytic residues. One A-loop conformer forms lasting interactions with the L16 segment, leading to reduced RMSF and greater compaction in the active site. By contrast, simulations of 0P-ERK2 reveal excursions of A-loop residues away from the C-lobe, leading to greater active site mobility. Thus, the A-loop in ERK2 switches between distinct conformations that reflect coupling with the active site, possibly via the L16 segment. Crystal packing interactions suggest that lattice contacts with the A-loop may restrain its structural variation in X-ray structures of ERK2. The novel conformational states identified by MD expand our understanding of ERK2 regulation, by linking the activated state of the kinase to reduced dynamics and greater compaction surrounding the catalytic site.

Proton transfer function of carbonic anhydrase: Insights from QM/MM simulations.

Recent QM/MM analyses of proton transfer function of human carbonic anhydrase II (CAII) are briefly reviewed. The topics include a preliminary analysis of nuclear quadrupole coupling constant calculations for the zinc ion and more detailed analyses of microscopic pK(a) of the zinc-bound water and free energy profile for the proton transfer. From a methodological perspective, our results emphasize that performing sufficient sampling is essential to the calculation of all these quantities, which reflects the well solvated nature of CAII active site. From a mechanistic perspective, our analyses highlight the importance of electrostatics in shaping the energetics and kinetics of proton transfer in CAII for its function. We argue that once the pK(a) for the zinc-bound water is modulated to be in the proper range (approximately 7.0), proton transfer through a relatively well solvated cavity towards/from the protein surface (His64) does not require any major acceleration. Therefore, although structural details like the length of the water wire between the donor and acceptor groups still may make a non-negligible contribution, our computational results and the framework of analysis suggest that the significance of such "fine-tuning" is likely secondary to the modulation of pK(a) of the zinc-bound water. We encourage further experimental analysis with mutation of (charged) residues not in the immediate neighborhood of the zinc ion to quantitatively test this electrostatics based framework; in particular, Phi analysis based on these mutations may shed further light into the relative importance of the classical Grotthus mechanism and the "proton hole" pathway that we have proposed recently for CAII.

Corrections to standard state in combustion calorimetry: an update and a web-based tool.

Combustion calorimetry is the predominant method for determination of enthalpies of formation for organic compounds. Both initial and final states of the calorimeter deviate significantly from the standard conditions. Correction of the obtained results to the standard state must be applied as accurately as possible to determine the combustion energy with an acceptable uncertainty, which is typically a few hundredths of a percent. The correction procedures in their current form were introduced in 1956 with simplifications to allow application in a pre-computer era. In this work, the procedures have been updated with respect to both the equations and reference values. The most reliable data sources are identified, and the updated algorithm is presented in the form of a Web-based tool available through the NIST TRC Web site.

Evaluating elastic network models of crystalline biological molecules with temperature factors, correlated motions, and diffuse x-ray scattering.

In this study, the variance-covariance matrix of protein motions is used to compare several elastic network models within the theoretical framework of x-ray scattering from crystals. A set of 33 ultra-high resolution structures is used to characterize the average scaling behavior of the vibrational density of states and make comparisons between experimental and theoretical temperature factors. Detailed investigations of the vibrational density of states, correlations, and predicted diffuse x-ray scatter are carried out for crystalline Staphylococcal nuclease; correlations and diffuse x-ray scatter are also compared to predictions from the translation, libration, screw model and a liquid-like dynamics model. We show that elastic network models developed to best predict temperature factors without regard for the crystal environment have relatively strong long-range interactions that yield very short-ranged atom-atom correlations. Further, we find that the low-frequency modes dominate the variance-covariance matrix only for those models with a physically reasonable vibrational density of states, and the fraction of modes required to converge the correlations is higher than that typically used for elastic network model studies. The practical implications are explored using computed diffuse x-ray scatter, which can be measured experimentally.

Application of elastic network models to proteins in the crystalline state.

Normal mode analysis using elastic network models has grown popular for probing the low-frequency collective dynamics of proteins and other biomolecular assemblies. In most previous studies, these models were validated by comparing calculated atomic fluctuations for isolated proteins with experimental temperature factors determined in the crystalline state, although there were also hints that including crystal contacts in the calculations has an impact on the comparison. In this study, a set of 83 ultra-high resolution crystal structures with experimentally determined anisotropic displacement parameters is used to evaluate several C(alpha)-based elastic network models that either ignore or treat the crystal environment in different ways; the latter include using periodic boundary conditions defined with respect to the asymmetric unit or the primitive unit cell as well as using the Born-von Kármán boundary condition that accounts for lattice vibrations. For all elastic network models, treating the crystal environment leads to better agreement with experimental anisotropic displacement parameters with the Born-von Kármán boundary condition giving the best agreement. Atomic correlations over the entire protein are clearly affected by the presence of the crystal contacts and fairly sensitive to the way that the crystal environment is treated. These observations highlight the importance of properly treating the protein system in an environment consistent with experiment when either evaluating approximate protein models or using approximate dynamic models in structural refinement application types. Finally, investigation of the scaling behaviors of the cumulative density of states and the heat capacity indicates that there are still gaps between simplified elastic models and all-atom models for proteins.

Interactions of the osmolyte glycine betaine with molecular surfaces in water: thermodynamics, structural interpretation, and prediction of m-values.

Noncovalent self-assembly of biopolymers is driven by molecular interactions between functional groups on complementary biopolymer surfaces, replacing interactions with water. Since individually these interactions are comparable in strength to interactions with water, they have been difficult to quantify. Solutes (osmolytes, denaturants) exert often large effects on these self-assembly interactions, determined in sign and magnitude by how well the solute competes with water to interact with the relevant biopolymer surfaces. Here, an osmometric method and a water-accessible surface area (ASA) analysis are developed to quantify and interpret the interactions of the remarkable osmolyte glycine betaine (GB) with molecular surfaces in water. We find that GB, lacking hydrogen bond donors, is unable to compete with water to interact with anionic and amide oxygens; this explains its effectiveness as an osmolyte in the Escherichia coli cytoplasm. GB competes effectively with water to interact with amide and cationic nitrogens (hydrogen bonding) and especially with aromatic hydrocarbon (cation-pi). The large stabilizing effect of GB on lac repressor-lac operator binding is predicted quantitatively from ASA information and shown to result largely from dehydration of anionic DNA phosphate oxygens in the protein-DNA interface. The incorporation of these results into theoretical and computational analyses will likely improve the ability to accurately model intra- and interprotein interactions. Additionally, these results pave the way for development of solutes as kinetic/mechanistic and thermodynamic probes of conformational changes and formation/disruption of molecular interfaces that occur in the steps of biomolecular self-assembly processes.

Validation of thermophysical data for scientific and engineering applications.

High quality thermophysical property data are essential to many scientific and engineering applications. These data are produced at a high rate and are affected by a range of experimental and reporting error sources that often exceed stated uncertainties. As a result, critical evaluation is required to establish the limits of reliability in a quantified way. The present work describes reporting recommendations and property data validation methods developed and applied at the Thermodynamics Research Center at NIST through the use of the ThermoData Engine (TDE; SRD 103a/b) software. Examples are provided with an emphasis on various consistency checks, which may include the use of equations of state (EOS).

Ensemble docking to difficult targets in early-stage drug discovery: Methodology and application to fibroblast growth factor 23.

Ensemble docking is now commonly used in early-stage in silico drug discovery and can be used to attack difficult problems such as finding lead compounds which can disrupt protein-protein interactions. We give an example of this methodology here, as applied to fibroblast growth factor 23 (FGF23), a protein hormone that is responsible for regulating phosphate homeostasis. The first small-molecule antagonists of FGF23 were recently discovered by combining ensemble docking with extensive experimental target validation data (Science Signaling, 9, 2016, ra113). Here, we provide a detailed account of how ensemble-based high-throughput virtual screening was used to identify the antagonist compounds discovered in reference (Science Signaling, 9, 2016, ra113). Moreover, we perform further calculations, redocking those antagonist compounds identified in reference (Science Signaling, 9, 2016, ra113) that performed well on drug-likeness filters, to predict possible binding regions. These predicted binding modes are rescored with the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) approach to calculate the most likely binding site. Our findings suggest that the antagonist compounds antagonize FGF23 through the disruption of protein-protein interactions between FGF23 and fibroblast growth factor receptor (FGFR).

Development of effective quantum mechanical/molecular mechanical (QM/MM) methods for complex biological processes.

Motivated by the long-term goal of understanding vectorial biological processes such as proton transport (PT) in biomolecular ion pumps, a number of developments were made to establish combined quantum mechanical/molecular mechanical (QM/MM) methods suitable for studying chemical reactions involving significant charge separation in the condensed phase. These developments were summarized and discussed with representative problems. Specifically, free energy perturbation and boundary potential methods for treating long-range electrostatics were implemented to test the robustness of QM/MM results for protein systems. It was shown that consistent models with sufficient sampling were able to produce quantitatively satisfactory results, such as pKa for titritable groups in the interior of T4-lysozyme, while an inconsistent treatment of electrostatics or lack of sufficient sampling may produce incorrect results. Modifications were made to an approximate density functional theory (SCC-DFTB) to improve the description of proton affinity and hydrogen-bonding, which are crucial for the treatment of PT in polar systems. Test calculations on water autoionization showed clearly that both improvements are necessary for quantitatively reliable results. Finally, the newly established SCC-DFTB/MM-GSBP protocol was used to explore mechanistic issues in carbonic anhydrase (CA). Preliminary results suggest that PT in CA occurs mainly through short water wires containing two water molecules in a thermally activated fashion. Although longer water wires occur with similar frequencies, PT along those pathways, on average, has substantially higher barriers, a result not expected based on previous studies. The fluctuations of water molecules peripheral to the water wire were found to make a larger impact on the PT energetics compared to polar protein residues in the active site, which are largely pre-organized and therefore have less tendency to reorganize during the reaction.

A computationally identified compound antagonizes excess FGF-23 signaling in renal tubules and a mouse model of hypophosphatemia.

Fibroblast growth factor-23 (FGF-23) interacts with a binary receptor complex composed of α-Klotho (α-KL) and FGF receptors (FGFRs) to regulate phosphate and vitamin D metabolism in the kidney. Excess FGF-23 production, which causes hypophosphatemia, is genetically inherited or occurs with chronic kidney disease. Among other symptoms, hypophosphatemia causes vitamin D deficiency and the bone-softening disorder rickets. Current therapeutics that target the receptor complex have limited utility clinically. Using a computationally driven, structure-based, ensemble docking and virtual high-throughput screening approach, we identified four novel compounds predicted to selectively inhibit FGF-23-induced activation of the FGFR/α-KL complex. Additional modeling and functional analysis found that Zinc13407541 bound to FGF-23 and disrupted its interaction with the FGFR1/α-KL complex; experiments in a heterologous cell expression system showed that Zinc13407541 selectivity inhibited α-KL-dependent FGF-23 signaling. Zinc13407541 also inhibited FGF-23 signaling in isolated renal tubules ex vivo and partially reversed the hypophosphatemic effects of excess FGF-23 in a mouse model. These chemical probes provide a platform to develop lead compounds to treat disorders caused by excess FGF-23.

pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols.

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.

Importance of van der Waals Interactions in QM/MM Simulations.

The importance of accurately treating van der Waals interactions between the quantum mechanical (QM) and molecular mechanical (MM) atoms in hybrid QM/MM simulations has been investigated systematically. First, a set of van der Waals (vdW) parameters was optimized for an approximate density functional method, the self-consistent charge-tight binding density functional (SCC-DFTB) approach, based on small hydrogen-bonding clusters. The sensitivity of condensed phase observables to the SCC-DFTB vdW parameters was then quantitatively investigated by SCC-DFTB/MM simulations of several model systems using the optimized set and two sets of extreme vdW parameters selected from the CHARMM22 forcefield. The model systems include a model FAD molecule in solution and a solvated enediolate, and the properties studied include the radial distribution functions of water molecules around the solute (model FAD and enediolate), the reduction potential of the model FAD and the potential of mean force for an intramolecular proton transfer in the enediolate. Although there are noticeable differences between parameter sets for gas-phase clusters and solvent structures around the solute, thermodynamic quantities in the condensed phase (e.g., reduction potential and potential of mean force) were found to be less sensitive to the numerical values of vdW parameters. The differences between SCC-DFTB/MM results with the three vdW parameter sets for SCC-DFTB atoms were explained in terms of the effects of the parameter set on solvation. The current study has made it clear that efforts in improving the reliability of QM/MM methods for energetical properties in the condensed phase should focus on components other than van der Waals interactions between QM and MM atoms.