The conversion of proteins between internal and cartesian coordinates is a limiting step in many pipelines, such as molecular dynamics simulations and machine learning models. This conversion is typically carried out by sequential or parallel applications of the Natural extension of Reference Frame (NeRF) algorithm. This work proposes a massively parallel NeRF implementation which, depending on the polymer length, achieves speedups between 400 and 1200× over the previous state-of-the-art. It accomplishes this by dividing the conversion into three main phases: parallel composition of the monomer backbone, assembly of backbone subunits, and parallel elongation of sidechains; and by batching these computations into a minimal number of efficient matrix operations. Special emphasis is placed on reusability and ease of use. We open source the code (available at https://github.com/EleutherAI/mp_nerf) and provide a corresponding python package.

In order to quantitatively predict nano- as well as other particle-size distributions, one needs to have both a mathematical model and estimates of the parameters that appear in these models. Here, we show how one can use Bayesian inversion to obtain statistical estimates for the parameters that appear in recently derived mechanism-enabled population balance models (ME-PBM) of nanoparticle growth. The Bayesian approach addresses the question of "how well do we know our parameters, along with their uncertainties?." The results reveal that Bayesian inversion statistical analysis on an example, prototype Ir(0)n nanoparticle formation system allows one to estimate not just the most likely rate constants and other parameter values, but also their SDs, confidence intervals, and other statistical information. Moreover, knowing the reliability of the mechanistic model's parameters in turn helps inform one about the reliability of the proposed mechanism, as well as the reliability of its predictions. The paper can also be seen as a tutorial with the additional goal of achieving a "Gold Standard" Bayesian inversion ME-PBM benchmark that others can use as a control to check their own use of this methodology for other systems of interest throughout nature. Overall, the results provide strong support for the hypothesis that there is substantial value in using a Bayesian inversion methodology for parameter estimation in particle formation systems.

For larger molecules, the computational demands of configuration selective vibrational configuration interaction theory (cs-VCI) are usually dominated by the configuration selection process, which commonly is based on second order vibrational Møller-Plesset perturbation (VMP2) theory. Here we present two techniques, which lead to substantial accelerations of such calculations while retaining the desired high accuracy of the final results. The first one introduces the concept of configuration classes, which allows for a highly efficient exploitation of the analogs of the Slater-Condon rules in vibrational structure calculations with large correlation spaces. The second approach uses a VMP2 like vector for augmenting the targeted vibrational wavefunction within the selection of configurations and thus avoids any intermediate diagonalization steps. The underlying theory is outlined and benchmark calculations are provided for highly correlated vibrational states of several molecules.

In the current work, the authors analyzed and compared the chemical behavior of some (anti)aromatic compounds. The species selected are benzene and cyclobutadiene as the aromatic and antiaromatic classical examples, respectively. Next, the anion Al42- , which is the first all-metal molecule catalogued as aromatic and its non-metallic isoelectronic analog, B42- were also analyzed. The antiaromatic clusters Al44- and B44- were studied in form of lithium salts. And the end, the non-planar B12 boron cluster and its isoelectronic analogs ( B11Be- and B11C+ ) were considered for being analyzed under the same criterium. The analysis was realized employing the following descriptors: molecular electrostatic potential and local hyper-softness to get insights concerning local reactivity when facing reagents leading to ionic or covalent interactions. The results show that all the molecules analyzed presented some specific variations in their respective local reactivity despite being labeled as aromatic compounds. This analysis provides a notion that the local reactivity is more based on the nature of the atoms and the molecular geometry than the aromatic character by itself.

Despite advances in the field, hemoincompatibility remains a critical issue for hemodialysis (HD) as interactions between various human blood constituents and the polymeric structure of HD membranes results in complications such as activation of immune system cascades. Adding hydrophilic polymer structures to the membranes is one modification approach that can decrease the extent of protein adsorption. This study conducted molecular dynamics (MD) simulations to understand the interactions between three human serum proteins (fibrinogen [FB], human serum albumin, and transferrin) and common HD membranes in untreated and modified forms. Poly(aryl ether sulfone) (PAES) and cellulose triacetate were used as the common dialyzer polymers, and membrane modifications were performed with 2-hydroxymethyl methacrylate (HEMA) and poly (2-methoxyethyl acrylate) (PMEA), using polydopamine-assisted co-deposition. The MD simulations were used as the framework for binding energy simulations, and molecular docking simulations were also performed to conduct molecular-level investigations between the two modifying polymers (HEMA and PMEA) and FB. Each of the three proteins acted differently with the membranes due to their unique nature and surface chemistry. The simulations show PMEA binds less intensively to FB with a higher number of hydrogen bonds, which reflects PMEA's superior performance compared to HEMA. The simulations suggest PAES membranes could be used in modified forms for blood-contact applications as they reflect the lowest binding energy to blood proteins.

We determine the proton affinity (PA) and gas-phase basicity (GB) of amino acid α-alanine at a chemically accurate level by performing explicitly-correlated CCSD(T)-F12b/aug-cc-pVDZ geometry optimizations and normal mode vibrational frequency calculations as well as CCSD(T)-F12b/aug-cc-pVTZ energy computations at the possible neutral and protonated geometries. Temperature effects at 298.15 K considering translational, rotational, and vibrational enthalpy and entropy corrections are obtained via standard statistical mechanics utilizing the molecular geometries and the harmonic vibrational energy levels. Both the amino nitrogen (N) and the carbonyl oxygen (O) atoms are proven to be potential protonation sites and a systematic conformational search reveals 3 N- and 9 O-protonated conformers in the 0.00-7.88 and 25.43-30.43 kcal/mol energy ranges at 0 K, respectively. The final computed PA and GB values at (0)298.15 K in case of N-protonation are (214.47)216.80 and 207.07 kcal/mol, respectively, whereas the corresponding values for O-protonation are (189.04)190.63 and 182.31 kcal/mol. The results of the benchmark high-level coupled-cluster computations are utilized to assess the accuracy of several lower-level cost-effective methods such as MP2 and density functional theory with various functionals (SOGGA11-X, M06-2X, PBE0, B3LYP, M06, TPSS).

The engineered three-helix bundle, UVF, is thermostabilized entropically due to heightened, native-state dynamics. However, it is unclear whether this thermostabilization strategy is observed in natural proteins from thermophiles. We performed all-atom, explicit solvent molecular dynamics simulations of two three-helix bundles from thermophilic H. butylicus (2lvsN and 2lvsC) and compared their dynamics to a mesophilic three-helix bundle, the Engrailed homeodomain (EnHD). Like UVF, 2lvsC had heightened native dynamics, which it maintained without unfolding at 100°C. Shortening and rigidification of loops in 2lvsN and 2lvsC and increased surface hydrogen bonds in 2lvsN were observed, as is common in thermophilic proteins. A buried disulfide and salt bridge in 2lvsN and 2lvsC, respectively, provided some stabilization, and addition of a homologous disulfide bond in EnHD slowed unfolding. The transferability and commonality of stabilization strategies among members of the three-helix bundle fold suggest that these strategies may be general and deployable in designing thermostable proteins.

AutoDock Vina (Vina) achieved a very high docking-success rate, p^ , but give a rather low correlation coefficient, R , for binding affinity with respect to experiments. This low correlation can be an obstacle for ranking of ligand-binding affinity, which is the main objective of docking simulations. In this context, we evaluated the dependence of Vina R coefficient upon its empirical parameters. R is affected more by changing the gauss2 and rotation than other terms. The docking-success rate p^ is sensitive to the alterations of the gauss1, gauss2, repulsion, and hydrogen bond parameters. Based on our benchmarks, the parameter set1 has been suggested to be the most optimal. The testing study over 800 complexes indicated that the modified Vina provided higher correlation with experiment Rset1 = 0.556 ± 0.025 compared with RDefault = 0.493 ± 0.028 obtained by the original Vina and RVina 1.2 = 0.503 ± 0.029 by Vina version 1.2. Besides, the modified Vina can be also applied more widely, giving R ≥ 0.500 for 32/48 targets, compared with the default package, giving R ≥ 0.500 for 31/48 targets. In addition, validation calculations for 1036 complexes obtained from version 2019 of PDBbind refined structures showed that the set1 of parameters gave higher correlation coefficient ( Rset1 = 0.617 ± 0.017 ) than the default package ( RDefault = 0.543 ± 0.020 ) and Vina version 1.2 ( RVina 1.2 = 0.540 ± 0.020 ). The version of Vina with set1 of parameters can be downloaded at https://github.com/sontungngo/mvina. The outcomes would enhance the ranking of ligand-binding affinity using Autodock Vina.

Molecular dynamics simulation is important in the computational study of the biomolecules. In this paper, we upgrade our previous FSATOOL to version 2.0. It is no longer a plugin as before. Besides the existed enhanced sampling and Markov state model analysis module, FSATOOL 2.0 has three new features now. First, it contains a molecular dynamics simulation engine on both CPU and GPU device. The engine works with an embedded enhanced sampling module. Second, it can do the free energy calculation by various practical methods, including the weighted histogram analysis method and Gaussian mixture model. Third, it has many subroutines to process the trajectory data, such as principal component analysis, time-structure based independent component analysis, contact analysis, and Φ-value analysis. Most importantly, all these calculations are integrated into one package. The trajectory data format is compatible with all the modules. With a proper input parameter file, users can do the molecular dynamics simulation and data analysis work by only a few simplified commands. The capabilities and theoretical backgrounds of FSATOOL 2.0 are introduced in the paper.

Thanks to its advantages, 19 F-NMR is an increasingly popular technique for the structural characterization of F-containing molecules, among which polymers, materials, fluorophores, pharmaceuticals, and so forth. However, the computational calculation of the 19 F-NMR chemical shifts, both for prediction and interpretation of experimental spectra, remains a challenge. In this work a density functional theory (DFT) based protocol for the calculation of the chemical shifts is established within the framework of the gauge-independent atomic orbital method, upon verifying the performance of Hartree-Fock and 60 DFT functionals coupled with seven different basis sets. The benchmark is conducted using two sets of molecules, namely one used for testing methods and another used for probing; the former set consists of 134 molecules, the latter 50, yet both of them with F in different chemical environments. Following Bally-Rablen-Tantillo strategy, the scaling parameters and other statistical quantities were computed for each method upon least squares linear regression between experimental and computed chemical shifts. The designed computational workflow is computationally inexpensive and represents a significant improvement with respect to the current state of the art.

The effect of a static electric (E)-field and an unchirped and chirped laser pulse field on the cycl[3.3.3]azine molecule was investigated using next-generation quantum theory of atoms in molecules (NG-QTAIM). Despite the magnitude of the E-field of the laser pulses being an order of magnitude lower than for the static E-field, the variation of the energy gap between the lowest lying singlet (S1 ) and triplet (T1 ) excited states was orders of magnitude greater for the laser pulse than for the static E-field. Insights into the response of the electronic structure were captured by NG-QTAIM, where differences in the inverted singlet-triplet gap due to the laser pulses were significant larger compared to those induced by the static E-field. The response of the S1 and T1 excited states, as determined by NG-QTAIM, switched discontinuously between weak and strong chemical character for the static E-field. In contrast, the response to the laser pulses, determined by NG-QTAIM, is to induce a continuous range of chemical character, indicating the unique ability of the laser pulses to induce polarization effects in the form of "mixed" bond types. Our analysis demonstrates that NG-QTAIM is a useful tool for understanding the response to laser irradiation of the lowest-lying singlet S1 and triplet T1 excited states of emitters exhibiting thermally activated delayed fluorescence. The chirped laser pulse led to more frequent instances of the desired outcome of an inverted singlet-triplet gap than the unchirped pulse, indicating its usefulness as a tool to design more efficient organic light-emitting diode devices.

When different proposals exist (or can reasonably be formulated) for the size of the unit cell (in terms of number of atoms) and space group of crystalline compounds, a strategy for exploring with simulation methods the various cases and for investigating their relative stability must be defined. The optimization schemes of periodic quantum mechanical codes work in fact at fixed space group and number of atoms per unit cell, so that only the fractional coordinates of the atoms and the lattice parameters are optimized. A strategy is here presented, based on four standard tools, used synergistically and in sequence: (1) the optimization of inner coordinates and unit cell parameters; (2) the calculation of the vibrational frequencies not only at Γ , but also at a set of k →  points (in the example presented here they are eight, generated by a shrinking factor 2), looking for possible negative wavenumbers. The latter correspond to maxima, rather than minima, along the coordinate described by the corresponding normal mode; (3) the exploration of the total energy along the mode with negative wavenumber, looking for the minimum of the curve; (4) the identification of the new space group corresponding to the reduced symmetry resulting from the previous step. The strategy is illustrated with reference to the KMnF3 perovskite compound, for which many space groups are proposed in the literature, ranging from cubic (Pm3¯m) to tetragonal (P4/bm or I4/cm) and orthorhombic (Pnma and Cmcm) down to monoclinic (P21 /m). The corresponding primitive cells contain 5, 10, and 20 atoms in the various cases, and the point symmetry reduces from 48 to 4 operators. In nature, the KMnF3 phase transitions also include the magnetic phases. For simplicity, here we limit the analysis to the ones that take place between ferromagnetic phases, as they are sufficiently rich for illustrating the proposed strategy. As the total energy differences involved can be as small as, say, 10-50 μHartree, a high numerical accuracy at each one of the steps mentioned above is required. The present calculations, performed with the CRYSTAL code, by using an all electron basis set and the Hartree-Fock and B3LYP functionals, document such an accuracy. The energy difference between the tetragonal I4/cm and cubic Pm3¯m phases is 225 μHartree, with a volume reduction of 0.58 Å3 ; the differences between the orthorhombic and tetragonal phases are an order of magnitude smaller, being 23 μHartree and 0.06 Å3 for total energy and cell volume, respectively.

Stable equilibrium compounds containing a planar antiaromatic cyclooctatetraene (COT) ring are promising candidates for organic electronic devices such as organic semiconductor transistors. The planarization of COT by incorporation into rigid planar π-systems, as well as by oxidation or reduction has attracted considerable attention in recent years. Using dispersion-corrected density functional theory calculations, we explore an alternative approach of planarizing COT derivatives by adsorption onto graphene. We show that strong π-π stacking interactions between graphene and COT derivatives induce a planar structure with an antiaromatic central COT ring. In addition to being reversible, this strategy provides a novel approach for planarizing COT without the need for incorporation into a rigid structure, atomic substitution, oxidation, or reduction.

Biochemical simuflation and analysis play a significant role in systems biology research. Numerous software tools have been developed to serve this area. Using these tools for completing tasks, for example, stochastic simulation, parameter fitting and optimization, usually requires sufficient computational power to make the duration of completion acceptable. COPASI is one of the most powerful tools for quantitative simulation and analysis targeted at biological systems. It supports systems biology markup language and covers multiple categories of tasks. This work develops an open source package ParaCopasi for parallel COPASI tasks and investigates its performance regarding accelerations. ParaCopasi can be installed on platforms equipped with multicore CPU to exploit the cores, scaling from desktop computers to large scale high-performance computing clusters. More cores bring more performance. The results show that the parallel efficiency has a positive correlation with the total workload. The parallel efficiency reaches a level of at least 95% for both homogeneous and heterogenous tasks when computational workload is adequate. An example is illustrated by applicating this package in parameter estimation to calibrate a biochemical kinetics model.

Docking studies play a critical role in the current workflow of drug discovery. However, limitations may often arise through factors including inadequate ligand sampling, a lack of protein flexibility, scoring function inadequacies (e.g., due to metals, co-factors, etc.), and difficulty in retaining explicit water molecules. Herein, we present a novel CHARMM-based induced fit docking (CIFDock) workflow that can circumvent these limitations by employing all-atom force fields coupled to enhanced sampling molecular dynamics procedures. Self-guided Langevin dynamics simulations are used to effectively sample relevant ligand conformations, side chain orientations, crystal water positions, and active site residue motion. Protein flexibility is further enhanced by dynamic sampling of side chain orientations using an expandable rotamer library. Steps in the procedure consisting of fixing individual components (e.g., the ligand) while sampling the other components (e.g., the residues in the active site of the protein) allow for the complex to adapt to conformational changes. Ultimately, all components of the complex-the protein, ligand, and waters-are sampled simultaneously and unrestrained with SGLD to capture any induced fit effects. This modular flexible docking procedure is automated using CHARMM scripting, interfaced with SLURM array processing, and parallelized to use the desired number of processors. We validated the CIFDock procedure by performing cross-docking studies using a data set comprised of 21 pharmaceutically relevant proteins. Five variants of the CHARMM-based SWISSDOCK scoring functions were created to quantify the results of the final generated poses. Results obtained were comparable to, or in some cases improved upon, commercial docking program data.

The accuracy of correlation energy recovered by coupled cluster single-, double-, and perturbative triple-excitations, CCSD(T), has led to the method being considered the gold standard of computational chemistry. The application of CCSD(T) has been limited to medium-sized molecular systems due to its steep scaling with molecular size. The recent development of alternative domain-based local pair natural orbital coupled-cluster method, DLPNO-CCSD(T), has significantly broadened the range of chemical systems to which CCSD(T) level calculations can be applied. Condensed systems such as ionic liquids (ILs) have a large contribution from London dispersion forces of up to 150 kJ mol-1 in large-scale clusters. Ionic liquids show appreciable charge transfer effects that result in the increased valence orbital delocalization over the entire ionic network, raising the question whether the application of methods based on localized orbitals is reliable for these semi-Coulombic materials. Here the performance of DLPNO-CCSD(T) is validated for the prediction of correlation interaction energies of two data sets incorporating single-ion pairs of protic and aprotic ILs. DLPNO-CCSD(T) produced results within chemical accuracy with tight parameter settings and a non-iterative treatment of triple excitations. To achieve spectroscopic accuracy of 1 kJ mol-1 , especially for hydrogen-bonded ILs and those containing halides, the DLPNO settings had to be increased by two orders of magnitude and include the iterative treatment of triple excitations, resulting in a 2.5-fold increase in computational cost. Two new sets of parameters are put forward to produce the performance of DLPNO-CCSD(T) within chemical and spectroscopic accuracy.

In this article we use MP2 and CCSD(T) calculations for the A24 and S66 data sets to explore how midbond functions can be used to generate cost effective counterpoise corrected supramolecular interaction energies of noncovalent complexes. We use the A24 data set to show that the primary role of midbond functions is not to approach the complete basis set limit, but rather to ensure a balanced description of the molecules and the interaction region (unrelated to the basis set superposition error). The need for balance is a consequence of using atom centered basis sets. In the complete basis set limit, the error will disappear, but reaching the complete basis set limit is not feasible beyond small systems. For S66 we investigate the need for increasing the number of midbond centers. Results show that adding a second midbond center increases the accuracy, but the effect is secondary to changing the atom centered basis set. Further, by comparing calculations using the 3s3p2d1f1g midbond set with using aug-cc-pVDZ and aug-cc-pVTZ as midbond sets, we see that the requirements for the midbond set to be effective, is not just that it contains diffuse functions, but also that high angular momentum functions are included. By comparing two approaches for placing midbond centers we show that results are not particularly sensitive to placement as long as the placement is reasonable.

A method for averaging of NMR parameters by molecular dynamics (MD) has been derived from the method of statistical averaging in MD snapshots, benchmarked and applied to structurally dynamic interpretation of the 31 P NMR shift (δ31P ) in DNA phosphates. The method employs adiabatic dependence of an NMR parameter on selected geometric parameter(s) that is weighted by MD-calculated probability distribution(s) for the geometric parameter(s) (Ad-MD method). The usage of Ad-MD for polymers is computationally convenient when one pre-calculated structural dependence of an NMR parameter is employed for all chemically equivalent units differing only in dynamic behavior. The Ad-MD method is benchmarked against the statistical averaging method for δ31P in the model phosphates featuring distinctively different structures and dynamic behavior. The applicability of Ad-MD is illustrated by calculating 31 P NMR spectra in the Dickerson-Drew DNA dodecamer. δ31P was calculated with the B3LYP/IGLO-III/PCM(water) and the probability distributions for the torsion angles adjacent to the phosphorus atoms in the DNA phosphates were calculated using the OL15 force field.