In the linear-scaling divide-and-conquer (DC) electronic structure method, each subsystem is calculated together with the neighboring buffer region, the size of which affects the energy error introduced by the fragmentation in the DC method. The DC self-consistent field calculation utilizes a scheme to automatically determine the appropriate buffer region that is as compact as possible for reducing the computational time while maintaining acceptable accuracy (J. Comput. Chem. 2018, 39, 909). To extend the automatic determination scheme of the buffer region to the DC second-order Møller-Plesset perturbation (MP2) calculation, a scheme for estimating the subsystem MP2 correlation energy contribution from each atom in the buffer region is proposed. The estimation is based on the atomic orbital Laplace MP2 formalism. Based on this, an automatic buffer determination scheme for the DC-MP2 calculation is constructed and its performance for several types of systems is assessed.

Accurate numerical calculations of porosities and related properties are of importance when analyzing metal-organic frameworks (MOFs). We present porE, an open-source, general-purpose implementation to compute such properties and discuss all results regarding their sensitivity to numerical parameters. Our code combines the numerical efficiency of Fortran with the user-friendliness of Python. Three different approaches to calculate porosities are implemented in porE, and their advantages and drawbacks are discussed. In contrast to commonly used implementations, our approaches are entirely deterministic and do not require any stochastic averaging. In addition to the calculation of porosities, porE can calculate pore size distributions and offers the possibility to analyze pore windows. The underlying approaches are outlined, and pore windows are discussed concerning their impact on the analyzed porosities. Comparisons with reference values aim for a clear differentiation between void and accessible porosities, which we provide for a small benchmark set consisting of eight MOFs. In addition, our approaches are used for a bigger benchmark set containing 370 MOFs, where we determine linear relationships within our approaches as well as to reference values. We show how these relationships can be used to derive corrections to a give porosity approach, minimizing its mean error. As a highlight we show how complex workflows can be designed with a few lines of Python code using porE.

Dynamical properties are of great importance in determining the behavior of synthetic and natural molecules, but capturing them by computational methods is a nontrivial task. Very often the time scales of the relevant phenomena are far beyond the typical time windows accessible by classical Molecular Dynamics (MD) simulations, currently limited to the order of microseconds on standard laboratory workstations. On the other hand, biased and accelerated simulations allow for fast and thorough exploration of the molecular conformational space, but they lose the dynamic information. The problem of recovering dynamics from biased/accelerated simulations is a very active field of research, but no totally robust/reliable solutions have been given yet. In this paper it is shown how the Smoluchowski equation, in the framework of Diffusion Theory (DT), can be used to bridge this gap, and dynamical properties, in the form of time correlation functions (TCFs), can be extracted also from such kind of simulations. DT is first extended (EDT) to express the mobility tensors entering the Smoluchowski operator in terms of a recently introduced unified and regularized Rotne-Prager-Yamakawa approximation, [P. J. Zuk, E. Wajnryb, K. A. Mizerski, P. Szymczak, J. Fluid. Mech. 2014, 741, R5, 1-13] also involving mixed rotation-translation contributions, and rotation-rotation terms beside the classical translation-translation ones, so far used in DT. Then, the method is applied to recover the dynamics of a nontrivial example of a peptide in explicit water from the first 200 ns of a Replica Exchange Molecular Dynamics simulation, which is a popular computational method that destroys the long time dynamics. EDT dynamics were found to favorably compare against those coming from a standard MD simulation of the same system, requiring a time window of 30 μs to converge. This result shows that EDT is a tool of practical value to recover the long time dynamics of systems in diffusive regimes from biased/accelerated simulations, to be exploited in those cases when direct evaluation by standard MD is unfeasible.

We introduce iSPECTRON, a program that parses data from common quantum chemistry software (NWChem, OpenMolcas, Gaussian, Cobramm, etc.), produces the input files for the simulation of linear and nonlinear spectroscopy of molecules with the Spectron code, and analyzes the spectra with a broad range of tools. Vibronic spectra are expressed in term of the electronic eigenstates, obtained from quantum chemistry computations, and vibrational/bath effects are incorporated in the framework of the displaced harmonic oscillator model, where all required quantities are computed at the Franck-Condon point. The program capabilities are illustrated by simulating linear absorption, transient absorption and two dimensional electronic spectra of the pyrene molecule. Calculations at two levels of electronic structure theory, time-dependent density functional theory (with NWChem) and RASSCF/RASPT2 (with OpenMolcas) are presented and compared where possible. The iSPECTRON program is available online at https://github.com/ispectrongit/iSPECTRON/ and distributed open source under the terms of the Educational Community License version 2.0 (ECL 2.0).

Two proteins within the β-grasp superfamily, the B1-domain of protein G and the small archaeal modifier protein 1, were investigated to elucidate the key determinants of structural stability at the level of individual interactions. These symmetrical proteins both contain two β-hairpins which form a sheet flanked by a central α-helix. They were subjected to high temperature molecular dynamics simulations and the detailed behavior of each long-range interaction was characterized. The results revealed that in GB1 the most stable region was the C-terminal hairpin and in SAMP1 it was the opposite, the N-terminal hairpin. Experimental results for GB1 support this finding. In conclusion, it appears that the difference in the location and number of hydrophobic interactions dictate the differential stability which is accommodated due to structural symmetry of the β-grasp fold. Thus, the hairpins are interchangeable and in nature this lends itself to adaptability and flexibility.

A detailed computational exploration of the most relevant steps of iodido Pt(IV) complexes reduction and Pt(II) drugs mechanism of action and eventual deactivation is presented here inspired by the recent findings on iodido Pt(II) complexes and surprising re-evaluation of their cytotoxic activity. Pt(II) and Pt(IV) model systems are investigated and compared with cisplatin and its Pt(IV) derivative. Both monodeprotonated ascorbic acid and l-cysteine are used as reducing agents in the inner-sphere reduction mechanism of Pt(IV) complexes. Aquation mechanism of iodido Pt(II) complexes, interaction with guanine and sulfur containing compounds and reaction with the model protein hen egg white lysozyme are explored, due to a detected different behavior with respect to classical platinum drugs. The outcomes of such exploration allow to shed light on the role that the increased soft character together with bridging and leaving abilities of iodide over chloride could play in determining the cytotoxic profile of iodido Pt drugs.

Efficient design and screening of the novel molecules is a major challenge in drug and material design. This paper focuses on a multi-stage pipeline, in which several deep neural network models are combined to map discrete molecular representations into continuous vector space to later generate from it new molecular structures with desired properties. Here, the Attention-based Sequence-to-Sequence model is added to "spellcheck" and correct generated structures, while the oversampling in the continuous space allows generating candidate structures with desired distribution for properties and molecular descriptors, even for a small reference datasets. We further use computer simulation to validate the desired properties in the numerical experiment. With the focus on the drug design, such a pipeline allows generating novel structures with a control of Synthetic Accessibility Score and a series of metrics that assess the drug-likeliness. Our code is available at https://github.com/SoftServeInc/novel-molecule-generation.

Excited state properties are difficult to trace back to the common molecular orbital picture when the excited state wavefunction is a linear combination of two or more Slater determinants. Here, a theoretical methodology is introduced based on the algebraic diagrammatic construction scheme for the polarization propagator (ADC(n)) that allows to make this connection and to eventually derive structure-function relationships. The usefulness of this approach is demonstrated by an analysis of the transition dipole moments of the low-lying 1B3u and 2B3u states of anthracene and (1,4,5,8)-tetraazaanthracene.

A more efficient version of the Expanded Ensembles method for calculation of free energy in molecular-mechanical simulations is proposed. The method is based on the Horowitz L2MC approach to accelerate movement along the "alchemical" coordinate. It is possible to achieve the same efficiency of the algorithm both with the optimal number of "windows" and with a larger number of them compared to the original algorithm. Since the optimal number of windows is unknown a priory, the proposed algorithm is more robust than the traditional one. We can choose the number of windows in excess and do not worry about the loss of efficiency. We illustrate the method's efficiency with the computation of the hydration free energies of pyridine and water.

Molecular dynamics (MD) simulations that rely on force field methods has been widely used to explore the structure and function of RNAs. However, the current commonly used force fields are limited by the electrostatic description offered by atomic charge, dipole and at most quadrupole moments, failing to capture the anisotropic picture of electronic features. Actually, the distribution of electrons around atomic nuclei is not spherically symmetric but is geometry dependent. A multipolar electrostatic model based on high rank multipole moments is described in this work, which allows us to combine polarizability and anisotropy of electron density. RNA secondary structure was taken as a research system, and its substructures including stem, loops (hairpin loop, bulge loop, internal loop, and multi-branch loop), and pseudoknots (H-type and K-type) were investigated, respectively, as well as the hairpin. First, the atom-atom electrostatic properties derived from one chain of a duplex RNA 2MVY in our previous work (Ref. 58) were measured by the pilot RNA systems of hairpin, hairpin loop, stem, and H-type pseudoknot, respectively. The prediction results were not satisfactory. Consequently, to obtain a general set of electrostatic parameters for RNA force fields, the convergence behavior of the atom-atom electrostatic interactions in the pilot RNA systems was explored using high rank atomic multipole moments. The pilot RNA systems were cut into four types of different-sized molecular fragments, and the single nucleotide fragment and nucleotide-paired fragment proved to be the most reasonable systems for base-unpairing regions and base-pairing regions to investigate the convergence behavior of all types of atom-atom electrostatic interactions, respectively. Transferability of the electrostatic properties drawn from the pilot RNA systems to the corresponding test systems was also investigated. Furthermore, the convergence behavior of atomic electrostatic interactions in other substructures including bulge loop, internal loop, multi-branch loop, and K-type pseudoknot was expected to be modeled via the hairpin.

Minimum-energy path (MEP) calculations, such as those typified by the nudged elastic band method, require input of reactant and product molecular configurations at initialization. In the case of reactions involving more than one molecule, generating initial reactant and product configurations requires careful consideration of the relative position and orientations of the reactive molecules in order to ensure that the resulting MEP calculation proceeds without converging on an alternative reaction-path, and without requiring excessive numbers of optimization iterations; as such, this initial system set-up is most commonly performed "by hand," with an expert user arranging reactive molecules in space to ensure that the following MEP calculation runs smoothly. In this Article, we introduce a simple preconditioning scheme which replaces this labor-intensive, human-knowledge-based step with an automated deterministic computational scheme. In our approach, initial reactant and product configurations are generated such that steric hindrance between reactive molecules is minimized in the reactant and product configurations, while also simultaneously requiring minimal structural differences between the reactants and products. The method is demonstrated using a benchmark test-set of >3400 organic molecular reactions, where comparison of the reactant/product configurations generated using our approach compare very well to initial configurations which were generated on an ad hoc basis.

This study aims to apply artificial neural networks for the prediction of the lattice parameters of materials with stannite- and kesterite-type structure, and to compare the results of predictions with that obtained in the calculations exploiting the density functional theory. Crystallographic data for 49 compounds with stannite-type structure and for four compounds with the kesterite-type structure are found and, based on it, crystal structures are calculated using the density functional theory (DFT) method in a two-step relaxation procedure for all compounds. An multilayer Perceptron is constructed, which then is trained on gathered crystallographic data. Values predicted by a neural network (lattice parameters) are compared with experimental data and with results of DFT calculations. Moreover, a global optimization method (the Uspex code) is used to find potentially novel crystal structures for investigated chemical compositions. The results are discussed in the term of advantages and disadvantages of each method.

We have performed a hierarchical ab initio benchmark and DFT performance study of D2 Ch•••A- chalcogen bonds (Ch = S, Se; D, A = F, Cl). The ab initio benchmark study is based on a series of ZORA-relativistic quantum chemical methods [HF, MP2, CCSD, CCSD(T)], and all-electron relativistically contracted variants of Karlsruhe basis sets (ZORA-def2-SVP, ZORA-def2-TZVPP, ZORA-def2-QZVPP) with and without diffuse functions. The highest-level ZORA-CCSD(T)/ma-ZORA-def2-QZVPP counterpoise-corrected complexation energies (ΔECPC ) are converged within 1.1-3.4 kcal mol-1 and 1.5-3.1 kcal mol-1 with respect to the method and basis set, respectively. Next, we used the ZORA-CCSD(T)/ma-ZORA-def2-QZVPP (ΔECPC ) as reference data for analyzing the performance of 13 different ZORA-relativistic DFT approaches in combination with the Slater-type QZ4P basis set. We find that the three-best performing functionals are M06-2X, B3LYP, and M06, with mean absolute errors (MAE) of 4.1, 4.2, and 4.3 kcal mol-1 , respectively. The MAE for BLYP-D3(BJ) and PBE amount to 8.5 and 9.3 kcal mol-1 , respectively.

We report six new dynamically stable structures of SrTiO3 at various pressures ranging from 0 to 200 GPa. These structures were found by exploring the enthalpy surface with the Minima Hopping structure prediction method. The potential energy surface was generated by a machine learned potential, the charge equilibration via neural network technique (CENT), based on an extensive training data set of highly diverse SrTiO3 periodic and cluster structures. All our CENT structures were validated at the level of density functional theory. For our new structures, we performed phonon calculations and NVT molecular dynamics calculations to investigate their dynamical stability. Finally, X-ray diffraction patterns were simulated to help to identify our predicted structures in experiments.

We report the application of our fragment-based quantum chemistry model MIM (Molecules-In-Molecules) with electrostatic embedding. The method is termed "EE-MIM (Electrostatically Embedded Molecules-In-Molecules)" and accounts for the missing electrostatic interactions in the subsystems resulting from fragmentation. Point charges placed at the atomic positions are used to represent the interaction of each subsystem with the rest of the molecule with minimal increase in the computational cost. We have carefully calibrated this model on a range of different sizes of clusters containing up to 57 water molecules. The fragmentation methods have been applied with the goal of reproducing the unfragmented total energy at the MP2/6-311G(d,p) level. Comparative analysis has been carried out between MIM and EE-MIM to gauge the impact of electrostatic embedding. Performance of several different parameters such as the type of charge and levels of fragmentation are analyzed for the prediction of absolute energies. The use of background charges in subsystem calculations improves the performance of both one- and two-layer MIM while it is noticeably important in the case of one-layer MIM. Embedded charges for two-layer MIM are obtained from a full system calculation at the low-level. For one-layer MIM, in the absence of a full system calculation, two different types of embedded charges, namely, Geometry dependent (GD) and geometry independent (GI) charges, are used. A self-consistent procedure is employed to obtain GD charges. We have further tested our method on challenging charged systems with stronger intermolecular interactions, namely, protonated ammonia clusters (containing up to 30 ammonia molecules). The observations are similar to water clusters with improved performance using embedded charges. Overall, the performance of NPA charges as embedded charges is found to be the best.

A high-speed numerical potential delivering computational performance comparable with complex coarse-grained analytic potentials makes available models that have greater degrees of physical and chemical accuracy. This opens the possibility of increased accuracy in classical molecular dynamics simulations of anisotropic systems. In this work, we report the development of a high-speed lookup table (LUT) of four-dimensional gridded data, that uses cubic B-spline interpolations to derive off grid values and their associated partial derivatives that are located between the known grid data points. The accuracy of the coarse-grained numerical potential using a LUT from uniaxial Gay-Berne (GB) potential produced array of values is within a 3% and a 5% margin of error respectively for the interpolation of the uniaxial GB potential and its partial derivatives. The numerical potential model and partial derivatives speedup is made competitive with the analytical potential by exploiting graphics processing units on board functionality. The capability of the numerical potential is demonstrated by comparing minimizations of a box of 500 naphthalene molecules. The minimizations using a full atomistic (NAMD/CHARMM force field), a biaxial GB and a numerical potential from a LUT using data from the CHARMM pair potential was done. The numerical potential model is significantly more accurate in its approximation of the atomistic local minimum configuration than is the biaxial GB analytical potential function. This demonstrates that using a numerical potential founded on a direct lookup of the atomistic potential landscape significantly improves coarse grain (CG) modeling of complex molecules, possibly paving the way for accurate anisotropic system CG modeling.

In this study we aim to determine the origin of the electron density describing a CH···HC interaction in planar and twisted conformers of biphenyl. In order to achieve this, the fragment, atomic, localized, delocalized, intra- and inter-atomic (FALDI) decomposition scheme was utilized to decompose the density in the inter-nuclear region between the ortho-hydrogens in both conformers. Importantly, the structural integrity, hence also topological properties, were fully preserved as no 'artificial' partitioning of molecules was implemented. FALDI-based qualitative and quantitative analysis revealed that the majority of electron density arises from two, non-classical and non-local effects: strong overlap of ortho CH σ-bonds, and long-range electron delocalization between the phenyl rings and ortho carbons and hydrogens. These effects resulted in a delocalized electron channel, that is, a density bridge or a bond path in a QTAIM terminology, linking the H-atoms in the planar conformer. The same effects and phenomena are present in both conformers of biphenyl. We show that the CH···HC interaction is a molecular-wide event due to large and long-range electron delocalization, and caution against approaches that investigate CH···HC interactions without fully taking into account the remainder of the molecule.

Among the different noncovalent interactions, halogen bonds have captured wide attention in the last years. Their stability has been rationalized in electrostatic terms by appealing to the σ-hole concept, a charge-depleted region that is able to interact favorably with electron rich moieties. This interpretation has been questioned, and in this work a set of anionic halogen model systems are used to shed some light on this issue. We use the interacting quantum atoms method, which provides an orbital invariant energy decomposition in which pure electrostatic terms are well isolated, and we complement our insights with the analysis of electrostatic potentials (ESPs) as well as with traditional descriptors of charge accumulation like the Laplacian of the electron density. The total electrostatic interaction between the interacting species is surprisingly destabilizing in many of the systems examined, demonstrating that although σ-holes might be qualitatively helpful, much care has to be taken in ascribing the stability of these systems to electrostatics. It is clearly shown that electron delocalization is essential to understand the stability of the complexes. The evolution of atomic charges as the aggregates forms reveals a charge transfer picture in which the central, σ-hole bearing halogen acts as a mere spectator. These systems may then be not far from engaging in a classical 3c-4e interaction. Since the presence of a σ-hole as characterized by the ESP mapped on a suitable molecular envelope isosurface does not guarantee attractive electrostatic interactions, we encourage to employ a wider perspective that takes into account the full charge distribution.