First-principles data set of 45,892 isolated and cation-coordinated conformers of 20 proteinogenic amino acids

We present a structural data set of the 20 proteinogenic amino acids and their amino-methylated and acetylated (capped) dipeptides. Different protonation states of the backbone (uncharged and zwitterionic) were considered for the amino acids as well as varied side chain protonation states. Furthermore, we studied amino acids and dipeptides in complex with divalent cations (Ca2+, Ba2+, Sr2+, Cd2+, Pb2+, and Hg2+). The database covers the conformational hierarchies of 280 systems in a wide relative energy range of up to 4 eV (390 kJ/mol), summing up to a total of 45,892 stationary points on the respective potential-energy surfaces. All systems were calculated on equal first-principles footing, applying density-functional theory in the generalized gradient approximation corrected for long-range van der Waals interactions. We show good agreement to available experimental data for gas-phase ion affinities. Our curated data can be utilized, for example, for a wide comparison across chemical space of the building blocks of life, for the parametrization of protein force fields, and for the calculation of reference spectra for biophysical applications.


Background & Summary
Proteins are the machinery of life. We here present a first-principles study of the conformational preferences of their basic building blocks-specifically, as summarized in Fig. 1: 20 proteinogenic amino acids and dipeptides, with different possible protonation states, and the conformational space changes resulting from attaching six divalent cations, i.e., Ca 2+ , Ba 2+ , Sr 2+ , Cd 2+ , Pb 2+ , and Hg 2+ . In past studies, a wide range of different approximate electronic structure methods has been applied to some of these proteinogenic amino acids-see, for example, references  . These studies have deepened our understanding of the conformational basics of individual building blocks, but a systematic comparison of properties of the different building blocks is complicated when relying on data from different sources. On the one hand this is due to the molecular models that may differ in protonation states and backbone capping. On the other, the simulations can differ in several ways: • Different sampling strategies or methods to generate conformers may have been used.
Search-dependent settings, like energy cut-offs, can also have a significant impact on the results. • The levels of theory that have been applied range from semi-empirical to Hartree-Fock (HF) to density-functional theory (DFT) up to coupled-cluster calculations  . • Numerical settings, e.g., basis sets, can differ substantially and might lead to different results.
A further point that limits a quantitative comparison is the accessibility of the data from different studies. Energies, for example, often have to be extracted from table footnotes and/or the structural data is not always accessible in the Supplementary Information of the respective articles, sometimes even only accessible as figures in the manuscript. The data set presented here overcomes such limitations by covering a comprehensive segment of chemical space exhaustively, using a large scale computational effort. This study treats 20 proteinogenic amino acids, their dipeptides and their interactions with the divalent cations Ca 2+ , Ba 2+ , Sr 2+ , Cd 2+ , Pb 2+ , and Hg 2+ (see Fig. 1 for an overview) on the same theoretical footing. The importance of peptide cation interactions may be highlighted by the fact that about 40% of all proteins bind cations [60][61][62] . Especially Ca 2+ is important in a multitude of functions, ranging,  for example, from blood clotting 63 to cell signaling to bone growth 64 . Such calcium mediated functions can be disturbed by the presence of alternative divalent heavy metal cations like Pb 2+ , Cd 2+ , and Hg 2+ (refs 62,65,66).
The conformations and total energies of each molecular system are calculated from first principles in the framework of density-functional theory (DFT) 67,68 using the PBE generalized-gradient exchange-correlation functional 69 . Energies are corrected for van der Waals interactions using the Tkatchenko-Scheffler formalism 70 . In this formalism, pairwise C 6 [n]/r 6 terms are computed and summed up for all pairs of atoms. r is the interatomic distance, a cut-off for short interatomic distances is applied, and C 6 [n] coefficients are obtained from the self-consistent electron density. The combined approach is referred to as 'PBE+vdW' throughout this work. This level of theory is robust for potential-energy surface (PES) sampling of peptide systems [71][72][73][74][75][76][77][78] . The curated data is provided as basis for comparative studies across chemical space to reveal conformational trends and energetic preferences. It can, for example, further be used for force-field development, theoretical studies at higher levels of theory, and as a starting point for theoretical calculations of spectra for biophysical applications.

Molecular models
This study covers a total of 280 molecular systems (summarized in Fig. 1). The number is the product of the following chemical degrees of freedom that were considered in our study: 20 proteinogenic amino acids. In case of (de)protonatable side chains, all protomers (different protonations states) were considered as well.
2 different backbone types, either free termini (considered in uncharged or zwitterionic form) or capped (N-terminally acetylated or C-terminally amino-methylated).
7 reflecting that the respective amino acid or dipeptide was considered either in isolation or with one of six different cation additions: Ca 2+ , Ba 2+ , Sr 2+ , Cd 2+ , Pb 2+ , or Hg 2+ .

Conformational search and energy functions
For the initial scan of the PES, the empirical force field OPLS-AA 79 was employed, followed by DFT-PBE+vdW relaxations of the energy minima identified in the force field. The identified set of structures was then subjected to a further first-principles refinement step, ab initio replica-exchange molecular dynamics (REMD). An overview of the procedure is given in Fig. 2 and the steps are described in more detail below.
Force-field based (OPLS-AA) 79 global conformational searches (Step 1) were performed for all dipeptides and amino acids (i) without a coordinating cation and (ii) with Ca 2+ . These searches employed a basin hopping search strategy 80,81 as implemented in the tool 'scan', distributed with the TINKER molecular simulation package 82,83 . We here use an in-house parallelized version of the TINKER scan utility that was first used in reference 74 . In this search strategy, input structures for relaxations are generated by projecting along normal modes starting from a local minimum. The number of search directions from a local minimum was set to 20. Conformers were accepted within a relative energy window of 100 kcal/mol and if they differ in energy from already found minima by at least 10 − 4 kcal/mol. The search terminates when the relaxations of input structures do not result in new minima.
After that, PBE+vdW relaxations (Step 2) were performed with the program FHI-aims [84][85][86] . FHI-aims employs numeric atom-centered orbital basis sets as described in reference 84 to discretize the Kohn-Sham orbitals. Different levels of computational defaults are available, distinguished by choice of the basis set, integration grids, and the order of the multipole expansion of the electrostatic (Hartree) potential of the electron density. For the chemical elements relevant to this work, 'light' settings include the so-called tier1 basis sets and were used for initial relaxations. 'Tight' settings include the larger tier2 basis sets and ensure converged conformational energy differences at a level of few meV (ref. 84). Unless noted otherwise, all energies discussed here are results of PBE+vdW calculations with a tier2 basis and 'tight' settings. Relativistic effects were taken into account by the so-called atomic zero-order regular approximation (atomic ZORA) 87,88 as described in reference 84 . Previous comparisons to high-level quantum chemistry benchmark calculations at the coupled-cluster level, CCSD(T), demonstrated the reliability of this approach for polyalanine systems 72,76 , alanine, phenylalanine, and glycine containing tripeptides 76 , and alanine dipeptides with Li + (ref. 73). Further benchmarks at the MP2 level of theory are reported below in the section Technical Validation.
The refinement (Step 3) by ab initio REMD 89,90 is intended to alleviate the potential effects of conformational energy landscape differences between the force field and the DFT method. In REMD, multiple molecular dynamics trajectories of the same system are independently initialized and run in a range of different temperatures. Based on a Metropolis criterion, configurations are swapped between trajectories of neighboring temperatures. Thus, the simulations can overcome barriers and provide an enhanced conformational sampling in comparison to classical molecular dynamics (MD) 90,91 . The simulations were carried out employing a script-based REMD scheme that is provided with FHI-aims and that was first used in reference 92  potential-energy surface as shown by Sindhikara et al. 93 . The velocity-rescaling approach by Bussi et al. 94 was used to sample the canonical distribution. Starting geometries for the replicas were taken from the lowest energy conformers resulting from the PBE+vdW relaxations in Step 2. REMD parameters for the individual systems, i.e. the number of replicas, acceptance rates for exchanges between replicas, the frequency for exchange attempts, and the temperature range, are summarized in Supplementary Table S1 of the Supplementary Information. Conformations were extracted from the REMD trajectories every 10th step, i.e. every 10 fs of simulation time. In order to generate a set of representative conformers, these structures were clustered using a k-means clustering algorithm 95 with a cluster radius of 0.3 Å as provided by the MMSTB package 96 . The resulting arithmetic-mean structures from each cluster were then relaxed using PBE+vdW with 'light' computational settings. The obtained conformers were again clustered and cluster representatives were relaxed with PBE+vdW ('tight' computational settings) to obtain the final conformation hierarchies. The refinement step by REMD is essential, as shown in Fig. 3, which separately identifies the number of distinct conformers found in Step 2 and, subsequently, the number of additional conformers found in Step 3.
After step 2, a total of 17,381 stationary points was found for the amino acids and dipeptides in isolation and in complex with Ca 2+ . The refinement procedure in Step 3 increases this number to a total of 21,259 structures. Initial structures for the Ba 2+ , Cd 2+ , Hg 2+ , Pb 2+ and Sr 2+ binding amino acid and dipeptide systems were then obtained by replacing the Ca 2+ cation in the amino acid and dipeptide structures binding a Ca 2+ cation. These structures were subsequently relaxed with PBE+vdW employing 'tight' computational settings and a tier-2 basis set. This procedure results in 24,633 further conformers with bound cations. Altogether, we thus provide information on 45,892 stationary points of the PBE+vdW PES for all systems studied in this work.
The numbers of conformers identified in the searches are also given in Supplementary Table S2 of the  Supplementary Information. Supplementary Tables S3 and S4 provide detailed accounts of how many structures were found for which amino acid/dipeptide in isolation or with attached cations.

Data Records
The curated data, consisting of the Cartesian coordinates of 45,892 stationary point geometries of the PBE +vdW PES (the main outcome of our work) and their potential energies computed at the 'tight'/tier-2 level of accuracy in the FHI-aims code, is provided as plain text files sorted in directories (see Fig. 4). The PBE+vdW total energies are included since they are an integral part of the construction of our geometry data sets. Importantly, the stationary point geometries could be used as starting points to refine the total energy accuracy by higher-level methods, e.g., those discussed in 'Technical Validation' below. The folder structure is hierarchic and straightforward. The naming scheme is explained in the following: control.in FHI-aims input file with technical parameters for the calculations. Please note that these files also include the exact specifications of the 'tight' numerical settings for all included elements. hierarchy_PBE+vdW_tier-2.dat in each final subfolder, contains three columns: number of the conformer, total energy (in eV, PBE+vdW, tier-2 basis set, 'tight' numerical settings, computed with FHI-aims version 031011), and relative energy (in eV, relative to the respective global minimum).
The curated data is publicly available from two sources: 1. A website dedicated to this data set has been set up (http://aminoaciddb.rz-berlin.mpg.de) and allows users to browse and download the data and to visualize molecular structures online.

Technical Validation
The conformational coverage for the amino acid alanine is validated by comparing to a recent study by Maul et al. 12 . In that reference 10    after global search and after REMD-refinement are shown in Fig. 5a. The results of our search (with the refinement step) are in good agreement with the data from reference 12 that is also shown in Fig. 5a. Structures are shown in Fig. 5b. Nine of the ten conformers identified by Maul et al. can be confirmed. The single conformer that is missing (highlighted by an X in Fig. 5a) is not a stationary point of the PBE+vdW potential energy surface. Conformers 14 and 15 are classified as saddle points by analysis of the vibrational modes.
In order to further quantify the reliability of the DFT-PBE+vdW level of theory for peptides, beyond earlier benchmark work 72,73,76 and especially with divalent cations, benchmark calculations were performed at the level of Møller-Plesset second-order perturbation theory (MP2) 99,100 using the electronic structure program package ORCA 101 . Single-point energy calculations were performed for all fixed stationary-point DFT-PBE+vdW geometries in our data base for the amino acids alanine (Ala) and phenylalanine (Phe) with neutral N and C termini in isolation as well as in complex with a Ca 2+ cation. Phe was selected to represent a 'difficult' example, i.e., the interaction of the cation with a larger aromatic side chain. The MP2 calculations did not include any frozen-core treatment (including semicore states is essential for Ca 2+ ) and were performed using Dunning's correlation-consistent polarized core-valence basis sets (cc-pCVnZ), with n = T/Q/5 denoting the triple-zeta, quadruple-zeta, and quintuple-zeta basis sets respectively 102 . The calculated SCF (Hartree-Fock) and MP2 correlation energies were then individually extrapolated to the complete basis set (CBS) limit as follows: For SCF energies, we used the extrapolation strategy proposed by Karton and Martin 103 : A, α, and the CBS-extrapolated energy E CBS SCF are parameters determined from a least-squares fitting algorithm applied individually for each conformer. For the MP2 correlation energies, an extrapolation scheme proposed by Truhlar 104 was applied: Again, B, β, and the CBS-extrapolated energy E CBS corr are parameters determined from a least-squares fitting algorithm as before. A detailed account of all numbers is given in the Supplementary Information (Supplementary Table S5). Mean absolute errors between the density-functional approximation (DFA) relative energies and the basis-set extrapolated MP2 relative energies were calculated as follows: where the index i runs over all N conformations of a given data set. ΔE i in principle denotes the energy difference between conformer i and the lowest-energy conformer of the set. The adjustable parameter c is used to shift the MP2 and DFA conformational hierarchies versus one another to obtain the lowest possible MAE, rendering the reported MAE value independent of the choice of any reference structure. Fig. 6a shows the corresponding obtained mean absolute errors (MAE) and maximal errors (max i |ΔE i DFA − ΔE i MP2 +c|) of different DFA calculations-performed with the FHI-aims code-with respect to benchmarks on the MP2 level obtained as described above. Within FHI-aims, the accuracy of integration grids and of the electrostatic potential was also verified by comparing 'tight' and 'really_tight' numerical settings, giving virtually identical results. The DFA level of theory of PBE+vdW shows a MAE well within chemical accuracy of 1 kcal=mol 43 meV for both structural sets of Ala and Phe; for Phe, the maximal error is 2 kcal=mol. We next applied a different long-range dispersion treatment, a recent  Table 1.   Fig. 6b compares the same set of DFAs to MP2 benchmark energy hierarchies. However, obtaining basis-set converged total energies of the same accuracy as for the isolated peptides by straightforward CBS extrapolation proved remarkably more difficult when Ca 2+ was involved. The reason is traced to the significant and slow-converging correlation contribution of the Ca 2+ semicore electrons, which leads to large and conformation dependent basis set superposition errors (BSSE). This problem was verified for MP2 calculations in the FHI-aims and ORCA codes, with several different basis set prescriptions 107 , and for CCSD(T) calculations. Standard DFAs, if sufficiently accurate, have a significant advantage in this respect since they are not subject to comparable numerical convergence problems. To yet arrive at reliable CBS-extrapolated MP2 conformational energy differences, we thus subjected the SCF and correlation energies of each Ca 2+ coordinated conformation to a counterpoise correction 108,109 to minimize the effect of BSSE on the Ca 2+ correlation energy contribution, prior to performing CBS extrapolation as described above. For the example of Ala+Ca 2+ and assuming rigid conformers, the BSSE is estimated as: E AlaþCa 2þ (Ala) represents the energy of Ala evaluated in the union of the basis sets on Ala and Ca 2+ , E Ala (Ala) represents the energy of Ala evaluated in the basis set on Ala, etc. The individual BSSE errors are then subtracted from the SCF and correlation energy respectively. Phe+Ca 2+ is treated equivalently. Complete numerical details are given in the Supplementary Information (Supplementary Table S6). Following this procedure, the MAE and maximal error values of various DFAs compared to MP2 are well within 1 kcal/mol for Ala+Ca 2+ . The PBE+vdW MAE for Phe+Ca 2+ amounts to just above~2 kcal/mol. The contributions from both the many-body dispersion and the hybrid PBE0 functional improve the MAE to just above 1 kcal/mol at to PBE0+MBD* level of theory. The maximum errors in the energy hierarchies between individual conformers are correspondingly larger. Overall, this assessment shows that our data base of conformer geometries constitutes, e.g., an excellent starting point for more exhaustive future benchmark work of new electronic structure methods for cation-peptide systems. For example, it would be very interesting to explore how F12 approaches, which address the correlation energy convergence problem explicitly, fare for a broad range of different Ca 2+ containing conformations of our peptides.
As a final validation, we compare the correlation of calculated gas-phase amino acid-Ca 2+ binding energies to the binding energy hierarchy found experimentally in a study by Ho et al. 110 . We calculate Figure 7. Comparison of the gas-phase binding energies of Ca 2+ to different amino acids calculated in this work (y axis) to the experimentally inferred hierarchy of gas-phase binding energies of Ca 2+ to different amino acids by Ho et al. 110 The amino acids are ordered along the x axis from the highest to lowest experimental Ca 2+ binding energy. Protonated and deprotonated Asp and Glu are not included among the experimental data and are here shown as predictions. E binding is high for deprotonated Asp and Glu since these forms of the amino acid would carry a negative charge.
www.nature.com/sdata/ SCIENTIFIC DATA | 3:160009 | DOI: 10.1038/sdata.2016.9 binding energy at the PES level as Energies E denote the PBE+vdW Born-Oppenheimer potential energies, including E amino acid of the lowest-energy conformers of the isolated amino acid and E complex of the same amino acid in complex with a Ca 2+ ion. Experimentally 110 , the gas-phase Ca 2+ affinities of 18 proteinogenic amino acids were determined by fragmenting Ca 2+ complexes with a combinatoric library of tripeptides at T ≈ 330 K, recording the mass spectrometric peak intensities of different fragmentation products. Quantitative average relative binding energies of Ca 2+ to different amino acids were thus inferred and can be compared to our findings, albeit with several important experiment-theory differences: (i) Entropy effects 73,75,111 should affect the specific complexes probed experimentally but cannot be included into the calculated numbers in the exact same way, (ii) structural differences (e.g., protonation, dimerized amino acids) between the fragments recorded in experiment and the amino acids covered here, (iii) experimental Ca 2+ affinities are not given for Asp and Glu because their gas-phase acidities, needed for data conversion, are not known. Fig. 7 compares the experimentally and theoretically inferred Ca 2+ binding affinities qualitatively. The x-axis reflects the experimental binding affinity energy hierarchy, arranging amino acids from left to right in order of decreasing binding affinity. The y axis shows calculated binding energies according to equation (5). Perfect correlation of the experimental and calculated hierarchies would imply a strictly monotonic decrease of calculated E binding values from left to right. This monotonic trend is not obeyed exactly; however, in view of the significant differences (i) and (ii) above, the qualitative agreement is quite striking. Normalized correlation coefficients between the experimental (1) and calculated (2) binding affinity data were calculated following the formula: with s 12 being the covariance of data sets and s i being the standard deviations of data sets i = 1,2. The result, correlation coefficients of r 12 = 0.979 or 0.909 for uncapped amino acids or dipeptides, respectively, also point to an overall remarkably good agreement. Finally, Fig. 7 also gives predicted E binding values for protonated (overall system charge +2) and deprotonated (overall system charge +1) Asp and Glu, reflecting the significant electrostatic attraction between cations and negatively charged (deprotonated) Asp and Glu side chains. The binding energy data sets are included as Supplementary Table S5.

Usage Notes
The present data contains stationary-point geometries (mainly minima, but also saddle points since no routine normal-mode analysis was performed) on the potential energy surface of the 20 proteinogenic amino acids and dipeptides, either isolated or in complex with a divalent cation (Ca 2+ , Ba 2+ , Sr 2+ , Cd 2+ , Pb 2+ , Hg 2+ ). The users of this dataset may find openbabel 112 (www.openbabel.org) to be a useful tool to convert FHI-aims and xyz files to other common file formats in chemistry.