On the accurate estimation of free energies using the jarzynski equality

The Jarzynski equality is one of the most widely celebrated and scrutinized nonequilibrium work theorems, relating free energy to the external work performed in nonequilibrium transitions. In practice, the required ensemble average of the Boltzmann weights of infinite nonequilibrium transitions is estimated as a finite sample average, resulting in the so-called Jarzynski estimator, . Alternatively, the second-order approximation of the Jarzynski equality, though seldom invoked, is exact for Gaussian distributions and gives rise to the Fluctuation-Dissipation estimator . Here we derive the parametric maximum-likelihood estimator (MLE) of the free energy considering unidirectional work distributions belonging to Gaussian or Gamma families, and compare this estimator to . We further consider bidirectional work distributions belonging to the same families, and compare the corresponding bidirectional to the Bennett acceptance ratio () estimator. We show that, for Gaussian unidirectional work distributions, is in fact the parametric MLE of the free energy, and as such, the most efficient estimator for this statistical family. We observe that and perform better than and , for unidirectional and bidirectional distributions, respectively. These results illustrate that the characterization of the underlying work distribution permits an optimal use of the Jarzynski equality. © 2018 Wiley Periodicals, Inc.     

A Density Functional Theory Study on Nonlinear Optical Properties of Double Cage Excess Electron Compounds: Theoretically Design M[Cu(Ag)@(NH3)n](M = Be,Mg and Ca; n = 1–3)

In this work, we investigated the nonlinear optical (NLO) properties of excess electron electride molecules of M[Cu(Ag)@(NH₃)_n](M = Be, Mg and Ca; n = 1–3) using density functional theory (DFT). This electride molecules consist of an alkaline-earth (Be, Mg and Ca) together with transition metal (Cu and Ag) doped in NH₃ cluster. The natural population analysis of charge and their highest occupied molecular orbital suggests that the M[Cu(Ag)@(NH₃)_n] compound has excess electron like alkaline-earth metal form double cage electrides molecules, which exhibit a large static first hyperpolarizability () (electron contribution part) and one of which owns a peak value of 216,938 (a.u.) for Be[Ag@(NH₃)₂] and vibrational harmonic first hyperpolarizability () (nuclear contribution part) values and the ratio of /, namely, η values from 0.02 for Be[Ag@(NH₃)] to 0.757 for Mg[Ag@(NH₃)₃]. The electron density contribution in different regions on values mainly come from alkaline-earth and transition metal atoms by first hyperpolarizability density analysis, and also explains the reason why values are positive and negative. Moreover, the frequency-dependent values β(−2ω,ω,ω) are also estimated to make a comparison with experimental measures. © 2018 Wiley Periodicals, Inc.     

Substitutional carbon defects in silicon: A quantum mechanical characterization through the infrared and Raman spectra

The infrared (IR) and Raman spectra of eight substitutional carbon defects in silicon are computed at the quantum mechanical level by using a periodic supercell approach based on hybrid functionals, an all electron Gaussian type basis set and the CRYSTAL code. The single substitutional C _s case and its combination with a vacancy (C _sV and C _sSiV) are considered first. The progressive saturation of the four bonds of a Si atom with C is then examined. The last set of defects consists of a chain of adjacent carbon atoms C, with i = 1–3. The simple substitutional case, C _s, is the common first member of the three sets. All these defects show important, very characteristic features in their IR spectrum. One or two C related peaks dominate the spectra: at 596 cm⁻¹ for C _s (and C _sSiV, the second neighbor vacancy is not shifting the C _s peak), at 705 and 716 cm⁻¹ for C _sV, at 537 cm⁻¹ for C and C (with additional peaks at 522, 655 and 689 for the latter only), at 607 and 624 cm⁻¹, 601 and 643 cm⁻¹, and 629 cm⁻¹ for SiC, SiC, and SiC, respectively. Comparison with experiment allows to attribute many observed peaks to one of the C substitutional defects. Observed peaks above 720 cm⁻¹ must be attributed to interstitial C or more complicated defects.     

A Computational Analysis of the Intrinsic Plasticity of Five-Coordinate Cu(II) Complexes and the Factors Leading to the Breakdown of the Orbital Directing Effect in Paddlewheel Secondary Building Units

Quantum chemical calculations on model copper paddlewheel (CPW) complexes of general formula [Cu₂(μ₂-O₂CR)₄L₂] establish two local coordination geometries at the metal centers depending on the balance between equatorial and axial ligand fields. When the equatorial field is stronger than the axial field (large ligand field asymmetry), dominates the stereochemical activity of the d⁹ shell resulting in a relatively rigid, “orbitally directed” planar or square pyramidal structure. However, if the axial field is significantly increased, or the equatorial field moderately weakened, a small ligand field asymmetry results and both and are involved in the stereochemical activity. This results in a “plastic,” distorted trigonal bipyramidal geometry where the former axial ligand moves into one of the original four equatorial positions. Linkers already used to synthesize zinc-dabco MOFs (dabco = 1,4-diazabicyclo[2.2.2]octane) are shown to generate plastic CPW secondary building unit analogs with potential implications for conferring breathing behavior for MOFs which would currently be assumed to be rigid. © 2019 Wiley Periodicals, Inc.     

Efficient singular-value decomposition of the coupled-cluster triple excitation amplitudes

We demonstrate a novel technique to obtain singular-value decomposition (SVD) of the coupled-cluster triple excitations amplitudes, . The presented method is based on the Golub-Kahan bidiagonalization strategy and does not require to be stored. The computational cost of the method is comparable to several coupled cluster singles and doubles (CCSD) iterations. Moreover, the number of singular vectors to be found can be predetermined by the user and only those singular vectors which correspond to the largest singular values are obtained at convergence. We show how the subspace of the most important singular vectors obtained from an approximate triple amplitudes tensor can be used to solve equations of the CC3 method. The new method is tested for a set of small and medium-sized molecular systems in basis sets ranging in quality from double- to quintuple-zeta. It is found that to reach the chemical accuracy (≈1 kJ/mol) in the total CC3 energies as little as 5 − 15% of SVD vectors are required. This corresponds to the compression of the amplitudes by a factor of about 0.0001 − 0.005 . Significant savings are obtained also in calculation of interaction energies or rotational barriers, as well as in bond-breaking processes. © 2019 Wiley Periodicals, Inc.     

Mono-silicon isoelectronic replacement in CAl4: van't hoff/le bel carbon or not?

In cluster studies, the isoelectronic replacement strategy has been successfully used to introduce new elements into a known structure while maintaining the desired topology. The well-known penta-atomic 18 valence electron (ve) species and its Al⁻/Si or Al/Si⁺ isoelectronically replaced clusters CAl₃Si⁻, CAl₂Si₂, , and , all possess the same anti-van't Hoff/Le Bel skeletons, that is, nontraditional planar tetracoordinate carbon (ptC) structure. In this article, however, we found that such isoelectronic replacement between Si and Al does not work for the 16ve-CAl₄ with the traditional van't Hoff/Le Bel tetrahedral carbon (thC) and its isoelectronic derivatives CAl₃X (X = Ga/In/Tl). At the level of CCSD(T)/def2-QZVP//B3LYP/def2-QZVP, none of the global minima of the 16ve mono-Si-containing clusters CAl₂SiX⁺ (X = Al/Ga/In/Tl) maintains thC as the parent CAl₄ does. Instead, X = Al/Ga globally favors an unusual ptC structure that has one long C─X distance yet with significant bond index value, and X = In/Tl prefers the planar tricoordinate carbon. The frustrated formation of thC in these clusters is ascribed to the CSi bonding that prefers a planar fashion. Inclusion of chloride ion would further stabilize the ptC of CAl₂SiAl⁺ and CAl₂SiGa⁺. The unexpectedly disclosed CAl₂SiAl⁺ and CAl₂SiGa⁺ represent the first type of 16ve-cationic ptCs with multiple bonds. © 2019 Wiley Periodicals, Inc.     

Fluxional bonds in quasi-planar and half-sandwich (M = K,Rb, and Cs)

Detailed molecular orbital and bonding analyses reveal the existence of both fluxional σ- and π-bonds in the global minima C_s ( 1 ) and C_s MB₁₈ ( 3 ) and transition states C_s ( 2 ) and C_s ( 4 ) of dianion and monoanions (M = K, Rb, and Cs). It is the fluxional bonds that facilitate the fluxional behaviors of the quasi-planar and half-sandwich which possess energy barriers smaller than the difference of the corresponding zero-point corrections. © 2019 Wiley Periodicals, Inc.     

Topological analysis of chemical bonding in the layered FePSe3 upon pressure-induced phase transitions

Two pressure-induced phase transitions have been theoretically studied in the layered iron phosphorus triselenide (FePSe₃ ). Topological analysis of chemical bonding in FePSe₃ has been performed based on the results of first-principles calculations within the periodic linear combination of atomic orbitals (LCAO) method with hybrid Hartree-Fock-DFT B3LYP functional. The first transition at about 6 GPa is accompanied by the symmetry change from to C2/m , whereas the semiconductor-to-metal transition (SMT) occurs at about 13 GPa leading to the symmetry change from C2/m to . We found that the collapse of the band gap at about 13 GPa occurs due to changes in the electronic structure of FePSe₃ induced by relative displacements of phosphorus or selenium atoms along the c-axis direction under pressure. The results of the topological analysis of the electron density and its Laplacian demonstrate that the pressure changes not only the interatomic distances but also the bond nature between the intralayer and interlayer phosphorus atoms. The interlayer P–P interactions are absent in two non-metallic FePSe₃ phases while after SMT the intralayer P–P interactions weaken and the interlayer P–P interactions appear.     

Accurate,robust, and reliable calculations of Poisson–Boltzmann binding energies

Poisson–Boltzmann (PB) model is one of the most popular implicit solvent models in biophysical modeling and computation. The ability of providing accurate and reliable PB estimation of electrostatic solvation free energy, , and binding free energy, , is important to computational biophysics and biochemistry. In this work, we investigate the grid dependence of our PB solver (MIBPB) with solvent excluded surfaces for estimating both electrostatic solvation free energies and electrostatic binding free energies. It is found that the relative absolute error of obtained at the grid spacing of 1.0 Å compared to at 0.2 Å averaged over 153 molecules is less than 0.2%. Our results indicate that the use of grid spacing 0.6 Å ensures accuracy and reliability in calculation. In fact, the grid spacing of 1.1 Å appears to deliver adequate accuracy for high throughput screening. © 2017 Wiley Periodicals, Inc.     

Ab initio structure and vibration-rotation dynamics of germylene,GeH2

Accurate structure and potential energy surface of germylene, GeH₂, in its ground electronic state ¹A₁ were determined from ab initio calculations using the coupled-cluster approach in conjunction with the correlation-consistent basis sets up to sextuple-zeta quality. The Born-Oppenheimer equilibrium structural parameters for the ¹A₁ state are estimated to be r_e(GeH) = 1.5793 Å and ∠_e(HGeH) = 91.19^∘. The term value T_e for the lowest excited electronic state ã³B₁ of GeH₂ is predicted to be 9140 cm^–1. The vibration-rotation energy levels for the ¹A₁ state of the ⁷⁴GeH₂, ⁷⁴GeD₂, ⁷²GeH₂, and ⁷⁰GeH₂ isotopologues were determined using a variational approach and compared with the experimental data. The role of the core-electron correlation, higher-order valence-electron correlation, scalar relativistic, spin-orbit, and adiabatic effects for prediction of the structure and vibration-rotation dynamics of the GeH₂ molecule is discussed. © 2019 Wiley Periodicals, Inc.     

Aspartate: An interesting model for analyzing dipole-ion and ion pair interactions through its oppositely charged amine and acid groups

Anionic species of aspartic acid, Asp⁻, having a zwitterionic backbone and a deprotonated side chain, appears to be a good example for analyzing dipole-ion and ion pair interactions. Density functional theory calculations were herein performed to investigate the low energy conformers of Asp⁻ embedded in a dielectric continuum modeling an aqueous environment, through a scan of the potential energy as a function of the side chain (χ₁, χ₂) torsion angles. The most energetically favorable conformers having g⁺g⁻ and g⁻g⁺ side chain orientations are found to be stabilized by charge-enhanced intramolecular H-bonding involving the positively charged () and the two negatively charged (COO⁻) groups. These conformers were further used to analyze Asp⁻ + nW clusters (W: water, n = 1 or 3), and Asp⁻/Asp⁻ pair formation. COO⁻ groups were found to be the most attractive sites for hosting a water molecule (binding energy: −6.0 ± 1.5 kcal/mol), compared to groups (binding energy: −4.7 ± 1.1 kcal/mol). Energy separation between g⁺g⁻ and g⁻g⁺ conformers increases upon explicit hydration. Asp⁻/Asp⁻ ion pairs, stabilized by the interaction between the group of a partner and the COO⁻ group of the other, shows a quite constant binding energy (−8.1 ± 0.2 kcal/mol), whatever the pair type, and the relative orientation of the two interacting partners. This study suggests a first step to achieve a more realistic image of intermolecular interactions in aqueous environment, especially upon increasing concentration. It can also be considered as a preliminary attempt to assess the interactions of the Lys⁺…Asp⁻/Glu⁻ ion pairs stabilizing intra- and interchain interactions in proteins.     

Update on phosphate and charged post‐translationally modified amino acid parameters in the GROMOS force field

In this study, we propose newly derived parameters for phosphate ions in the context of the GROMOS force field parameter sets. The non‐bonded parameters used up to now lead to a hydration free energy, which renders the dihydrogen phosphate ion too hydrophobic when compared to experimentally derived values, making a reparametrization of the phosphate moiety necessary. Phosphate species are of great importance in biomolecular simulations not only because of their crucial role in the backbone of nucleic acids but also as they represent one of the most important types of post‐translational modifications to protein side‐chains and are an integral part in many lipids. Our re‐parametrization of the free dihydrogen phosphate (H PO ) and three derivatives (methyl phosphate, dimethyl phosphate, and phenyl phosphate) leads, in conjunction with the previously updated charged side‐chains in the GROMOS parameter set 54A8, to new nucleic acid backbone parameters and a 54A8 version of the widely used GROMOS protein post‐translational modification parameter set. © 2017 Wiley Periodicals, Inc.     

The performance of density functional and wavefunction‐based methods for 2D and 3D structures of Au10

The transition from 2D to 3D structures in small gold clusters occurs around 10 atoms. Density functional theory predicts a planar structure for in contrast to recent second‐order Møller–Plesset perturbation theory calculations, which predict a 3D arrangement. The validity of the use of single‐reference second‐order Møller–Plesset theory for near metallic systems remains, however, questionable. On the other hand, it is less than clear how well density functional approximations perform for such clusters. We, therefore, decided to carry out quantum chemical calculations for using a variety of different density functionals as well as wavefunction‐based methods including coupled cluster theory to compare the different energetically low lying 2D and 3D cluster isomers. The results are perhaps not encouraging showing that most computational methods do not predict correctly the energetic sequence of isomers compared to coupled cluster theory. As perturbative triple corrections in the coupled cluster treatment change the order in cluster stability, the onset of 2D to 3D transition in these gold clusters remains elusive. As expected, second‐order Møller–Plesset theory is not suitable for correctly describing such systems.     

Mechanism of the O2(1Δg) generation from the Cl2/H2O2 basic aqueous solution explored by the combined ab initio calculation and nonadiabatic dynamics simulation

In the present work, mechanism of the O₂(¹Δ_g) generation from the reaction of the dissolved Cl₂ with H₂O₂ in basic aqueous solution has been explored by the combined ab initio calculation and nonadiabatic dynamics simulation, together with different solvent models. Three possible pathways have been determined for the O₂(¹Δ_g) generation, but two of them are sequentially downhill processes until formation of the OOCl⁻ complex with water, which are of high exothermic character. Once the complex is formed, singlet molecular oxygen is easily generated by its decomposition along the singlet-state pathway. However, triplet molecular oxygen of O₂() can be produced with considerable probability through nonadiabatic intersystem crossing in the ¹Δ_g/ intersection region. It has been found that the coupled solvent, heavy-atom, and nonadiabatic effects have an important influence on the quantum yield of the O₂(¹Δ_g) generation. © 2018 Wiley Periodicals, Inc.     

Development of a new parameter optimization scheme for a reactive force field based on a machine learning approach

Reactive molecular dynamics (MD) simulation is performed using a reactive force field (ReaxFF). To this end, we developed a new method to optimize the ReaxFF parameters based on a machine learning approach. This approach combines the k-nearest neighbor and random forest regressor algorithm to efficiently locate several possible ReaxFF parameter sets. As a pilot test of the developed approach, the optimized ReaxFF parameter set was applied to perform chemical vapor deposition (CVD) of an α-Al₂O₃ crystal. The crystal structure of α-Al₂O₃ was reasonably reproduced even at a relatively high temperature (2000 K). The reactive MD simulation suggests that the (110) surface grows faster than the (0001) surface, indicating that the developed parameter optimization technique could be used for understanding the chemical reaction in the CVD process. © 2019 Wiley Periodicals, Inc.     

Low‐energy structures of benzene clusters with a novel accurate potential surface

The benzene‐benzene (Bz‐Bz) interaction is present in several chemical systems and it is known to be crucial in understanding the specificity of important biological phenomena. In this work, we propose a novel Bz‐Bz analytical potential energy surface which is fine‐tuned on accurate ab initio calculations in order to improve its reliability. Once the Bz‐Bz interaction is modeled, an analytical function for the energy of the clusters may be obtained by summing up over all pair potentials. We apply an evolutionary algorithm (EA) to discover the lowest‐energy structures of clusters (for ), and the results are compared with previous global optimization studies where different potential functions were employed. Besides the global minimum, the EA also gives the structures of other low‐lying isomers ranked by the corresponding energy. Additional ab initio calculations are carried out for the low‐lying isomers of and clusters, and the global minimum is confirmed as the most stable structure for both sizes. Finally, a detailed analysis of the low‐energy isomers of the n = 13 and 19 magic‐number clusters is performed. The two lowest‐energy isomers show S₆ and C₃ symmetry, respectively, which is compatible with the experimental results available in the literature. The structures reported here are all non‐symmetric, showing two central Bz molecules surrounded by 12 nearest‐neighbor monomers in the case of the five lowest‐energy structures. © 2015 Wiley Periodicals, Inc.     

Software package SIMPRE—Revisited

This article elucidates the pitfalls identified in the software package SIMPRE recently developed by Baldoví et al. (J. Comput. Chem. 2013, 34, 1961) for modeling the spectroscopic and magnetic properties of single ion magnets as well as single‐molecule magnets. Analysis of the methodology used therein reveals that the crystal field parameters (CFPs), expressed nominally in the Stevens formalism, exhibit features characteristic for the CFPs expressed in the Wybourne notation. The resemblance of the two types of CFPs introduces a serious confusion that may lead to wrong comparisons of the CFPs taken from various sources. To clarify this confusion, the properties of the CFPs ( , ) associated with the Stevens operators ( X = S , J , or L ), which belong to the class of the tesseral‐tensor operators, are contrasted with those of the CFPs B_kq associated with the Wybourne operators , which belong to the class of the spherical‐tensor operators. Importantly, the confused properties of Stevens and Wybourne operators may bear on reliability of SIMPRE calculations. To consider this question independent calculations are carried out using the complete approach and compared with those of the restricted approach utilized earlier. It appears that the numerical results of the package SIMPRE are formally acceptable, however, the meaning of the CFPs must be properly reformulated. Several other conceptual problems arising from misinterpretations of the crucial notions and the CFP notations identified therein are also discussed and clarified. © 2014 Wiley Periodicals, Inc.     

Computational study of basis set and electron correlation effects on anapole magnetizabilities of chiral molecules

In the presence of a static, nonhomogeneous magnetic field, represented by the axial vector at the origin of the coordinate system and by the polar vector , assumed to be spatially uniform, the chiral molecules investigated in this paper carry an orbital electronic anapole, described by the polar vector . The electronic interaction energy of these molecules in nonordered media is a cross term, coupling and via , one third of the trace of the anapole magnetizability a_αβ tensor, that is, . Both and W^BC have opposite sign in the two enantiomeric forms, a fact quite remarkable from the conceptual point of view. The magnitude of predicted in the present computational investigation for five chiral molecules is very small and significantly biased by electron correlation contributions, estimated at the density functional level via three different functionals. © 2016 Wiley Periodicals, Inc.