The Poisson-Boltzmann equation is a widely used model to study electrostatics in molecular solvation. Its numerical solution using a boundary integral formulation requires a mesh on the molecular surface only, yielding accurate representations of the solute, which is usually a complicated geometry. Here, we utilize adjoint-based analyses to form two goal-oriented error estimates that allow us to determine the contribution of each discretization element (panel) to the numerical error in the solvation free energy. This information is useful to identify high-error panels to then refine them adaptively to find optimal surface meshes. We present results for spheres and real molecular geometries, and see that elements with large error tend to be in regions where there is a high electrostatic potential. We also find that even though both estimates predict different total errors, they have similar performance as part of an adaptive mesh refinement scheme. Our test cases suggest that the adaptive mesh refinement scheme is very effective, as we are able to reduce the error one order of magnitude by increasing the mesh size less than 20% and come out to be more efficient than uniform refinement when computing error estimations. This result sets the basis toward efficient automatic mesh refinement schemes that produce optimal meshes for solvation energy calculations.

Rare event methods applied to molecular simulations are growing in popularity, accessible and customizable software solutions have thus been developed and released. One of the most recent is PyRETIS, an open Python library for performing path sampling simulations. Here, we introduce PyVisA, a postprocessing package for path sampling simulations, which includes visualization and analysis tools for interpreting path sampling outputs. PyVisA integrates PyRETIS functionalities and aims to facilitate the determination of: (a) the correlation of the order parameter with other descriptors; (b) the presence of latent variables; and (c) intermediate meta-stable states. To illustrate some of the main PyVisA features, we investigate the proton transfer reaction in a protonated water trimer simulated via a simple polarizable model (Stillinger-David).

An algorithm to apply bond-angle constraints in molecular dynamics simulations of macromolecules or molecular liquids is presented. It uses Cartesian coordinates and determines the Lagrange multipliers required for maintaining the constraints iteratively. It constitutes an alternative to the use of only distance constraints (DCs) between particles to maintain a particular geometry. DCs are unsuitable to maintain particular, for example, linear or flat, geometries of molecules. The proposed algorithm can easily handle bond-length, bond-angle, and dihedral-angle constraints simultaneously, as when calculating a potential of mean force along a dihedral-angle degree of freedom.

QMCube (QM3 ) is a suite written in the Python programming language, initially focused on multiscale QM/MM simulations of biological systems, but open enough to address other kinds of problems. It allows the user to combine highly efficient QM and MM programs, providing unified access to a wide range of computational methods. The suite also supplies additional modules with extra functionalities. These modules facilitate common tasks such as performing the setup of the models or process the data generated during the simulations. The design of QM3 has been carried out considering the least number of external dependencies (only an algebra library, already included in the distribution), which makes it extremely portable. Also, the modular structure of the suite should help to expand and develop new computational methods.

Comment in J Comput Chem. 2021 Jul 15;42(19):1338-1340.

IOData is a free and open-source Python library for parsing, storing, and converting various file formats commonly used by quantum chemistry, molecular dynamics, and plane-wave density-functional-theory software programs. In addition, IOData supports a flexible framework for generating input files for various software packages. While designed and released for stand-alone use, its original purpose was to facilitate the interoperability of various modules in the HORTON and ChemTools software packages with external (third-party) molecular quantum chemistry and solid-state density-functional-theory packages. IOData is designed to be easy to use, maintain, and extend; this is why we wrote IOData in Python and adopted many principles of modern software development, including comprehensive documentation, extensive testing, continuous integration/delivery protocols, and package management. This article is the official release note of the IOData library.

The concept of effective one-electron potentials (EOPs) has proven to be extremely useful in efficient description of electronic structure of chemical systems, especially extended molecular aggregates such as interacting molecules in condensed phases. Here, a general method for EOP-based elimination of electron repulsion integrals is presented, that is tuned toward the fragment-based calculation methodologies such as the second generation of the effective fragment potentials (EFP2) method. Two general types of the EOP operator matrix elements are distinguished and treated either via the distributed multipole expansion or the extended density fitting (DF) schemes developed in this work. The EOP technique is then applied to reduce the high computational costs of the effective fragment charge-transfer (CT) terms being the bottleneck of EFP2 potentials. The alternative EOP-based CT energy model is proposed, derived within the framework of intermolecular perturbation theory with Hartree-Fock noninteracting reference wavefunctions, compatible with the original EFP2 formulation. It is found that the computational cost of the EFP2 total interaction energy calculation can be reduced by up to 38 times when using the EOP-based formulation of CT energy, as compared to the original EFP2 scheme, without compromising the accuracy for a wide range of weakly interacting neutral and ionic molecular fragments. The proposed model can thus be used routinely within the EFP2 framework.

Comment on J Comput Chem. 2020 May 15;41(13):1242-1251.

We are introducing Pep McConst-a software that employs a Monte-Carlo algorithm to construct 3D structures of polypeptide chains which could subsequently be studied as stand-alone macromolecules or complement the structure of known proteins. Using an approach to avoid steric clashes, Pep McConst allows to create multiple structures for a predefined primary sequence of amino acids. These structures could then effectively be used for further structural analysis and investigations. The article introduces the algorithm and describes its user-friendly approach that was made possible through the VIKING online platform. Finally, the manuscript provides several highlight examples where Pep McConst was used to predict the structure of the C-terminal of a known protein, generate a missing bit of already crystallized protein structures and simply generate short polypeptide chains.

We present, to the best of our knowledge, the first attempt to exploit the super-computer platform for quantum chemical density matrix renormalization group (QC-DMRG) calculations. We have developed the parallel scheme based on the in-house MPI global memory library, which combines operator and symmetry sector parallelisms, and tested its performance on three different molecules, all typical candidates for QC-DMRG calculations. In case of the largest calculation, which is the nitrogenase FeMo cofactor cluster with the active space comprising 113 electrons in 76 orbitals and bond dimension equal to 6000, our parallel approach scales up to approximately 2000 CPU cores.

The Ni-B complex 1BCF with a facilely accessible monophosphine (Pt Bu3 ) unit was theoretically designed, which was found to be more active than that with an ambiphilic ligand for hydrogenation of styrene. Substituting Pt Bu3 with a stronger electron donating ligand N-heterocyclic carbene largely improves the activity of the Ni-B complex.

A computational method to investigate the global conformational change of a protein is proposed by combining the linear response path following (LRPF) method and three-dimensional reference interaction site model (3D-RISM) theory, which is referred to as the LRPF/3D-RISM method. The proposed method makes it possible to efficiently simulate protein conformational changes caused by either solutions of varying concentrations or the presence of cosolvent species by taking advantage of the LRPF and 3D-RISM. The proposed method is applied to the urea-induced denaturation of ubiquitin. The LRPF/3D-RISM trajectories successfully simulate the early stage of the denaturation process within the simulation time of 300 ns, whereas no significant structural change is observed even in the 1 μs standard MD simulation. The obtained LRPF/3D-RISM trajectories reproduce the mechanism of the urea denaturation of ubiquitin reported in previous studies, and demonstrate the high efficiency of the method.

In this work the possibility of using the IR intensity of the stretching vibration νs of proton donor group for estimation of hydrogen bond strength was investigated. For a set of complexes with FH···X (X = F, N, O) hydrogen bonds in the wide range of energies (0.1-49.2 kcal/mol) vibrational frequencies νs and their intensities A were calculated (CCSD at complete basis set limit). The validity of the previously proposed linear proportionality between the intensification of the stretching vibration νs in IR spectra and hydrogen bond enthalpy -ΔH = 12.2 [Formula: see text] (A. V. Iogansen, Spectrochimica Acta A 1999) was examined. It is shown that for a range of similar hydrogen bond types with complexation energies ∆E <15 kcal/mol the ∆E( [Formula: see text] ) function remains similar to that proposed in the Iogansen's work, while upon strengthening this dependency becomes significantly nonlinear. We examined two other parameters ( [Formula: see text] and [Formula: see text] ) related to IR intensity as descriptors of hydrogen bond strength which are proportional to transition dipole moment matrix element and mass-independent dipole moment derivative. It was found that the dependency ∆E( [Formula: see text] ) stays linear in the whole studied range of complexation energies and it can be used for evaluation of ∆E from infrared spectral data with the accuracy about 2 kcal/mol. The mass-independent product [Formula: see text] is an appropriate descriptor for sets of complexes with various hydrogen bond types. Simple equations proposed in this work can be used for estimations of hydrogen bond strength in various systems, where experimental thermodynamic methods or direct calculations are difficult or even impossible.

We evaluate the correlation between binding energy (BE) and electron density ρ(r) at the bond critical point for 28 neutral hydrogen bonds, recently reported by Emamian and co-workers (J. Comput. Chem., 2019, 40, 2868). As an efficient tool, we use local stretching force constant kHBa derived from the local vibrational mode theory of Konkoli and Cremer. We compare the physical nature of BE versus kHBa , and provide an important explanation for cases with significant deviation in the BE- kHBa relation as well as in the BE-ρ(r) correlation. We also show that care has to be taken when different hydrogen bond strength measures are compared. The BE is a cumulative hydrogen bond strength measure while kHBa is a local measure of hydrogen bond strength covering different aspects of bonding. A simplified and unified description of hydrogen bonding is not always possible and needs an in-depth understanding of the systems involved.

To understand the effect of f-functions in predicting the right reaction mechanism for hypervalent iodine reagents, we adopt the Ahlrichs basis set family def2-SVP and def2-TZVP to revisit the potential energy surfaces of IBX-mediated oxidation and Togni I's isomerisation. Our results further prove that f-functions (in either Pople, Dunning, or Ahlrichs basis set series) are indispensable to predict the correct rate-determining step of hypervalent iodine reagents. The f-functions have a significant impact on the predicted reaction barriers for processes involving the IX (X = O, OH, CF3 , etc.) bond cleavage and formation, for example, in the reductive elimination step or the hypervalent twist step. We furthermore explore two hypervalent twist modes that account for the different influences of f-functions for IBX and Togni I. Our findings may be helpful for theoretical chemists to appropriately study the reaction mechanism of hypervalent iodine reagents.

Time-dependent (TD) density functional theory (DFT) and Franck-Condon Hertzberg-Teller (FCHT) calculations of various DFT functionals [B3LYP, CAM-B3LYP, ωB97XD, and optimally tuned (OT) long-range corrected (LC)-BLYP] were performed to examine how well DFT functionals can predict the experimental absorption and fluorescence spectra of a 12-carbon nanobelt (CNB). OT-LC-BLYP reproduced the experimental absorption spectrum well in terms of the peak position and intensity in the case of using a basis set with a diffuse function, such as 6-31+G(d,p) and 6-311+G(d,p), whereas B3LYP showed a red-shift in the peak positions and CAM-B3LYP and ωB97XD, which have a long-range HF exchange, showed blue shifts. Regarding the fluorescence spectrum calculations with FCHT using 6-311+G(d,p), the OT-LC-BLYP reproduced both the peak intensities and positions closest to the experimental spectrum. B3LYP, however, showed red-shifted peaks, and ωB97XD showed blue-shifted peaks. CAM-B3LYP provided less blue-shifted peaks, but the relative peak intensities mismatched the experimental ones. Furthermore, calculations of the absorption and vibrationally resolved fluorescence spectra of 16-CNB and 24-CNB using OT-LC-BLYP/6-311+G(d,p) showed absorption and fluorescence spectra close to the experimental spectra with high accuracy. Moreover, the application of a polarizable continuum model (dichloromethane) produced a red shift in the peak positions of the absorption spectrum with increasing intensity but an increase in the peak intensities of the fluorescence calculations without shifting the peak position.

We explore the possible route to approximate natural orbital occupation numbers-based diagnostic of differential multireference character of noncovalent energy differences by techniques based on many-body expansion. It turns out that two-body fragmentation of monomers may lead to a reasonable approximation of such a diagnostic in hydrogen-bonded complexes. The results are useful, for example, for assessment of the expected bias cancellation in energy differences of larger systems obtained by single-reference methods.

A detailed reaction mechanism of acetylene cyclotrimerization catalyzed by V(i PrNPMe2 )3 Fe-PMe3 (denote as CAT), a heterobimetallic complex featuring V-Fe triple bond, was computationally investigated using density functional theory. The calculated results show that the first acetylene firstly attaches to the V atom of CAT to get a four-membered ring structure through [2 + 2] cycloaddition reaction. For the second acetylene addition, there are two cyclotrimerization mechanisms, outer sphere mechanism and inner mechanism. The inner sphere reaction pathway is the main reaction pathway. By replacing the V with Nb and Ta, Fe with Ru and Os, a series of new catalysts are screened computationally. The calculated results show that, all of the nine heterobimetallic complexes show high activity at mild condition. The energy barrier of the rate determining step is related to the natural population analysis (NPA) charge of M' and the Wiberg bond index (WBI) of M-M' bond. The more negative NPA charge of M' and the smaller WBI of M-M' bond, the lower energy barrier is.

A local optimization algorithm for solving the Kohn-Sham equations is presented. It is based on a direct minimization of the energy functional under the equality constraints representing the Grassmann Manifold. The algorithm does not require an eigendecomposition, which may be advantageous in large-scale computations. It is optimized to reduce the number of Kohn-Sham matrix evaluations to one per iteration to be competitive with standard self-consistent field (SCF) approach accelerated by direct inversion of the iterative subspace (DIIS). Numerical experiments include a comparison of the algorithm with DIIS. A high reliability of the algorithm is observed in configurations where SCF iterations fail to converge or find a wrong solution corresponding to a stationary point different from the global minimum. The local optimization algorithm itself does not guarantee that the found minimum is global. However, a randomization of the initial approximation shows a convergence to the right minimum in the vast majority of cases.