Fractional Extended Diffusion Theory to capture anomalous relaxation from biased/accelerated molecular simulations.

Biased and accelerated molecular simulations (BAMS) are widely used tools to observe relevant molecular phenomena occurring on time scales inaccessible to standard molecular dynamics, but evaluation of the physical time scales involved in the processes is not directly possible from them. For this reason, the problem of recovering dynamics from such kinds of simulations is the object of very active research due to the relevant theoretical and practical implications of dynamics on the properties of both natural and synthetic molecular systems. In a recent paper [A. Rapallo et al., J. Comput. Chem. 42, 586-599 (2021)], it has been shown how the coupling of BAMS (which destroys the dynamics but allows to calculate average properties) with Extended Diffusion Theory (EDT) (which requires input appropriate equilibrium averages calculated over the BAMS trajectories) allows to effectively use the Smoluchowski equation to calculate the orientational time correlation function of the head-tail unit vector defined over a peptide in water solution. Orientational relaxation of this vector is the result of the coupling of internal molecular motions with overall molecular rotation, and it was very well described by correlation functions expressed in terms of weighted sums of suitable time-exponentially decaying functions, in agreement with a Brownian diffusive regime. However, situations occur where exponentially decaying functions are no longer appropriate to capture the actual dynamical behavior, which exhibits persistent long time correlations, compatible with the so called subdiffusive regimes. In this paper, a generalization of EDT will be given, exploiting a fractional Smoluchowski equation (FEDT) to capture the non-exponential character observed in the relaxation of intramolecular distances and molecular radius of gyration, whose dynamics depend on internal molecular motions only. The calculation methods, proper to EDT, are adapted to implement the generalization of the theory, and the resulting algorithm confirms FEDT as a tool of practical value in recovering dynamics from BAMS, to be used in general situations, involving both regular and anomalous diffusion regimes.

Extended diffusion theory: Recovering dynamics from biased/accelerated molecular simulations.

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.

Influence of the Ion Coordination Number on Cation Exchange Reactions with Copper Telluride Nanocrystals.

Cu2-xTe nanocubes were used as starting seeds to access metal telluride nanocrystals by cation exchanges at room temperature. The coordination number of the entering cations was found to play an important role in dictating the reaction pathways. The exchanges with tetrahedrally coordinated cations (i.e., with coordination number 4), such as Cd(2+) or Hg(2+), yielded monocrystalline CdTe or HgTe nanocrystals with Cu2-xTe/CdTe or Cu2-xTe/HgTe Janus-like heterostructures as intermediates. The formation of Janus-like architectures was attributed to the high diffusion rate of the relatively small tetrahedrally coordinated cations, which could rapidly diffuse in the Cu2-xTe NCs and nucleate the CdTe (or HgTe) phase in a preferred region of the host structure. Also, with both Cd(2+) and Hg(2+) ions the exchange led to wurtzite CdTe and HgTe phases rather than the more stable zinc-blende ones, indicating that the anion framework of the starting Cu2-xTe particles could be more easily deformed to match the anion framework of the metastable wurtzite structures. As hexagonal HgTe had never been reported to date, this represents another case of metastable new phases that can only be accessed by cation exchange. On the other hand, the exchanges involving octahedrally coordinated ions (i.e., with coordination number 6), such as Pb(2+) or Sn(2+), yielded rock-salt polycrystalline PbTe or SnTe nanocrystals with Cu2-xTe@PbTe or Cu2-xTe@SnTe core@shell architectures at the early stages of the exchange process. In this case, the octahedrally coordinated ions are probably too large to diffuse easily through the Cu2-xTe structure: their limited diffusion rate restricts their initial reaction to the surface of the nanocrystals, where cation exchange is initiated unselectively, leading to core@shell architectures. Interestingly, these heterostructures were found to be metastable as they evolved to stable Janus-like architectures if annealed at 200 °C under vacuum.

Peptide dynamics by molecular dynamics simulation and diffusion theory method with improved basis sets.

Improved basis sets for the study of polymer dynamics by means of the diffusion theory, and tests on a melt of cis-1,4-polyisoprene decamers, and a toluene solution of a 71-mer syndiotactic trans-1,2-polypentadiene were presented recently [R. Gaspari and A. Rapallo, J. Chem. Phys. 128, 244109 (2008)]. The proposed hybrid basis approach (HBA) combined two techniques, the long time sorting procedure and the maximum correlation approximation. The HBA takes advantage of the strength of these two techniques, and its basis sets proved to be very effective and computationally convenient in describing both local and global dynamics in cases of flexible synthetic polymers where the repeating unit is a unique type of monomer. The question then arises if the same efficacy continues when the HBA is applied to polymers of different monomers, variable local stiffness along the chain and with longer persistence length, which have different local and global dynamical properties against the above-mentioned systems. Important examples of this kind of molecular chains are the proteins, so that a fragment of the protein transthyretin is chosen as the system of the present study. This peptide corresponds to a sequence that is structured in ß-sheets of the protein and is located on the surface of the channel with thyroxin. The protein transthyretin forms amyloid fibrils in vivo, whereas the peptide fragment has been shown [C. P. Jaroniec, C. E. MacPhee, N. S. Astrof, C. M. Dobson, and R. G. Griffin, Proc. Natl. Acad. Sci. U.S.A. 99, 16748 (2002)] to form amyloid fibrils in vitro in extended ß-sheet conformations. For these reasons the latter is given considerable attention in the literature and studied also as an isolated fragment in water solution where both experimental and theoretical efforts have indicated the propensity of the system to form ß turns or α helices, but is otherwise predominantly unstructured. Differing from previous computational studies that employed implicit solvent, we performed in this work the classical molecular dynamics simulation on a realistic model solution with the peptide embedded in an explicit water environment, and calculated its dynamic properties both as an outcome of the simulations, and by the diffusion theory in reduced statistical-mechanical approach within HBA on the premise that the mode-coupling approach to the diffusion theory can give both the long-range and local dynamics starting from equilibrium averages which were obtained from detailed atomistic simulations.

Development of a semiempirical potential for simulations of thiol-gold interfaces. Application to thiol-protected gold nanoparticles.

A new semiempirical potential, based on density functional calculations and a bond-order Morse-like potential, is developed to simulate the adsorption behavior of thiolate molecules on non-planar gold surfaces, including relaxing effects, in a more realistic way. The potential functions include as variables the metal-molecule separation, vibrational frequencies, bending and torsion angles between several pairs of atom types and the coordination number of both the metal (Au) and thiolate groups. The potential was parameterized based on a set of density functional calculations of molecular adsorption in several surface sites (i.e. hollow, bridge, top, on-top Au adatom and the novel staple motif) for different crystalline facets, i.e. Au(111) and (100). Langevin dynamics simulations have been performed to study the capping effects of alkanethiolates molecules on Au nanoparticles in the range 1-4 nm. The simulation results reveal an enhancement of the coverage degree whilst the nanoparticles diameter decreases. A high surface disorder due to the strong S-Au bond was found, in very good agreement with very recent experimental findings [M. M. Mariscal, J. A. Olmos-Asar, C. Gutierrez-Wing, A. Mayoral and M. J. Yacaman, Phys. Chem. Chem. Phys., 2010, 12, 11785].

Clusters in strong polyelectrolyte solutions in the condensation theory approach.

The interaction free energy of parallel clusters of like-charged rod polyelectrolytes in solution is calculated in the framework of the extended condensation theory. For sufficiently high linear charge density of the polyelectrolyte, clustering takes place. The greater is the number of polyelectrolytes participating to the cluster, the smaller is the equilibrium interpolyelectrolyte distance, and the deeper is the corresponding free energy minimum. It is a counterintuitive organization due to the increasing of the counterion condensed charge and condensation volume, taking place as the polyelectyrolytes approach each other.

VARICELLA: a variable-cell direct space method for structure determination from powder diffraction data.

A direct space method for structure determination from powder diffraction data is proposed. Employing a hybrid Monte Carlo algorithm for generating the random conformations of a flexible molecular model, and sampling in a modified multicanonical statistical ensemble, it allows for variable cell parameters during an iterative search process. The acceptance-rejection criterion involves both a disagreement factor between the calculated and the experimental diffraction profiles and a modified crystal energy so that the space of tentative solutions can be widely explored while maintaining some physical meaningfulness of the proposals. Allowing the cell to be variable requires the zero shift to be treated as an optimizing parameter; this, in turn, requiring the disagreement factor to be based on the Fourier transform of the spectrum. The algorithm is presented in both a serial and a parallel version, the latter presenting several advantages, such as the possibility to probe different structures at a time while keeping them far from each other in the space defined by suitable order parameters. The method is built up and carefully tested by using, as a case study, a crystal of 3-ethyl 2,3-exo-disyndiotactic norbornene heptamer recently determined by single crystal x-ray diffraction techniques.

Formulation of improved basis sets for the study of polymer dynamics through diffusion theory methods.

In this work a new method is proposed for the choice of basis functions in diffusion theory (DT) calculations. This method, named hybrid basis approach (HBA), combines the two previously adopted long time sorting procedure (LTSP) and maximum correlation approximation (MCA) techniques; the first emphasizing contributions from the long time dynamics, the latter being based on the local correlations along the chain. In order to fulfill this task, the HBA procedure employs a first order basis set corresponding to a high order MCA one and generates upper order approximations according to LTSP. A test of the method is made first on a melt of cis-1,4-polyisoprene decamers where HBA and LTSP are compared in terms of efficiency. Both convergence properties and numerical stability are improved by the use of the HBA basis set whose performance is evaluated on local dynamics, by computing the correlation times of selected bond vectors along the chain, and on global ones, through the eigenvalues of the diffusion operator L. Further use of the DT with a HBA basis set has been made on a 71-mer of syndiotactic trans-1,2-polypentadiene in toluene solution, whose dynamical properties have been computed with a high order calculation and compared to the "numerical experiment" provided by the molecular dynamics (MD) simulation in explicit solvent. The necessary equilibrium averages have been obtained by a vacuum trajectory of the chain where solvent effects on conformational properties have been reproduced with a proper screening of the nonbonded interactions, corresponding to a definite value of the mean radius of gyration of the polymer in vacuum. Results show a very good agreement between DT calculations and the MD numerical experiment. This suggests a further use of DT methods with the necessary input quantities obtained by the only knowledge of some experimental values, i.e., the mean radius of gyration of the chain and the viscosity of the solution, and by a suitable vacuum trajectory, with great savings in computational time required. This offers a theoretical bridge between the experimental static and dynamical properties of polymers.

An algorithm for the uniform sampling of iso-energy surfaces and for the calculation of microcanonical averages.

In this article an algorithm is proposed to efficiently perform the uniform sampling of an iso-energy surface corresponding to a fixed potential energy U of a molecular system, and for calculating averages of certain quantities over microstates having this energy (microcanonical averages). The developed sampling technique is based upon the combination of a recently proposed method for performing constant potential energy molecular dynamics simulations [Rapallo, A. J Chem Phys 2004, 121, 4033] with well-established thermostatting techniques used in the framework of standard molecular dynamics simulations, such as the Andersen thermostat, and the Nose-Hoover chain thermostat. The proposed strategy leads to very accurate and drift-free potential energy conservation during the whole sampling process, and, very important, specially when dealing with high-dimensional or complicated potential functions, it does not require the calculation of the potential energy function hessian. The technique proved to be very reliable for sampling both low- and high-dimensional surfaces.

Global optimization of bimetallic cluster structures. II. Size-matched Ag-Pd, Ag-Au, and Pd-Pt systems.

Genetic algorithm global optimization of Ag-Pd, Ag-Au, and Pd-Pt clusters is performed. The 34- and 38-atom clusters are optimized for all compositions. The atom-atom interactions are modeled by a semiempirical potential. All three systems are characterized by a small size mismatch and a weak tendency of the larger atoms to segregate at the surface of the smaller ones. As a result, the global minimum structures exhibit a larger mixing than in Ag-Cu and Ag-Ni clusters. Polyicosahedral structures present generally favorable energetic configurations, even though they are less favorable than in the case of the size-mismatched systems. A comparison between all the systems studied here and in the previous paper (on size-mismatched systems) is presented.

Global optimization of bimetallic cluster structures. I. Size-mismatched Ag-Cu, Ag-Ni, and Au-Cu systems.

A genetic algorithm approach is applied to the optimization of the potential energy of a wide range of binary metallic nanoclusters, Ag-Cu, Ag-Ni, Au-Cu, Ag-Pd, Ag-Au, and Pd-Pt, modeled by a semiempirical potential. The aim of this work is to single out the driving forces that make different structural motifs the most favorable at different sizes and chemical compositions. Paper I is devoted to the analysis of size-mismatched systems, namely, Ag-Cu, Ag-Ni, and Au-Cu clusters. In Ag-Cu and Ag-Ni clusters, the large size mismatch and the tendency of Ag to segregate at the surface of Cu and Ni lead to the location of core-shell polyicosahedral minimum structures. Particularly stable polyicosahedral clusters are located at size N = 34 (at the composition with 27 Ag atoms) and N = 38 (at the composition with 32 and 30 Ag atoms). In Ag-Ni clusters, Ag32Ni13 is also shown to be a good energetic configuration. For Au-Cu clusters, these core-shell polyicosahedra are less common, because size mismatch is not reinforced by a strong tendency to segregation of Au at the surface of Cu, and Au atoms are not well accommodated upon the strained polyicosahedral surface.

Potential energy constrained molecular dynamics simulations.

A method for carrying out molecular dynamics simulations in which the potential energy U of the molecular system is constrained at its initial value is developed and thoroughly tested. The constraint is not introduced within the framework of the Lagrange multipliers technique, rather it is fulfilled in a natural way by carrying out the simulations in terms of suitable sets of delocalized coordinates. Such coordinates are defined by an appropriate tuning of the Baker, Kessi, and Delley internal delocalized nonredundant coordinates technique [J. Chem. Phys. 105, 192 (1996)]. The proposed method requires multiple evaluations of energy and gradients in each step of the molecular dynamics simulation, so that constant U simulations suffer some overhead compared to ordinary simulations. But the particular formulation of the delocalized coordinates and of the equations of motion greatly simplifies all the various steps required by the Baker's technique, thus allowing for the efficient implementation of the method itself. The technique is reliable and allows for very high accuracy in the potential energy conservation during the whole simulation. Moreover, it proved to be free of drift troubles which can occur when standard constraint methods are straightforwardly implemented without the application of appropriate correcting techniques.