Searching for globally optimal functional forms for interatomic potentials using genetic programming with parallel tempering

Molecular dynamics and other molecular simulation methods rely on a potential energy function, based only on the relative coordinates of the atomic nuclei. Such a function, called a force field, approximately represents the electronic structure interactions of a condensed matter system. Developing such approximate functions and fitting their parameters remains an arduous, time-consuming process, relying on expert physical intuition. To address this problem, a functional programming methodology was developed that may enable automated discovery of entirely new force-field functional forms, while simultaneously fitting parameter values. The method uses a combination of genetic programming, Metropolis Monte Carlo importance sampling and parallel tempering, to efficiently search a large space of candidate functional forms and parameters. The methodology was tested using a nontrivial problem with a well-defined globally optimal solution: a small set of atomic configurations was generated and the energy of each configuration was calculated using the Lennard-Jones pair potential. Starting with a population of random functions, our fully automated, massively parallel implementation of the method reproducibly discovered the original Lennard-Jones pair potential by searching for several hours on 100 processors, sampling only a minuscule portion of the total search space. This result indicates that, with further improvement, the method may be suitable for unsupervised development of more accurate force fields with completely new functional forms.     

Predicting stability limits for pure and doped dicationic noble gas clusters undergoing coulomb explosion: A parallel tempering based study

We have used a replica exchange Monte‐Carlo procedure, popularly known as Parallel Tempering, to study the problem of Coulomb explosion in homogeneous Ar and Xe dicationic clusters as well as mixed Ar–Xe dicationic clusters of varying sizes with different degrees of relative composition. All the clusters studied have two units of positive charges. The simulations reveal that in all the cases there is a cutoff size below which the clusters fragment. It is seen that for the case of pure Ar, the value is around 95 while that for Xe it is 55. For the mixed clusters with increasing Xe content, the cutoff limit for suppression of Coulomb explosion gradually decreases from 95 for a pure Ar to 55 for a pure Xe cluster. The hallmark of this study is this smooth progression. All the clusters are simulated using the reliable potential energy surface developed by Gay and Berne (Gay and Berne, Phys. Rev. Lett. 1982, 49, 194). For the hetero clusters, we have also discussed two different ways of charge distribution, that is one in which both positive charges are on two Xe atoms and the other where the two charges are at a Xe atom and at an Ar atom. The fragmentation patterns observed by us are such that single ionic ejections are the favored dissociating pattern. © 2017 Wiley Periodicals, Inc.     

Theory of slope-dependent disjoining pressure with application to Lennard-Jones liquid films

A liquid film of thickness h<100 nm is subject to additional intermolecular forces, which are collectively called disjoining pressure Pi. Since Pi dominates at small film thicknesses, it determines the stability and wettability of thin films. Current theory derived for uniform films gives Pi=Pi(h). This solution has been applied recently to non-uniform films and becomes unbounded near a contact line as h-->0. Consequently, many different effects have been considered to eliminate or circumvent this singularity. We present a mean-field theory of Pi that depends on the slope h(x) as well as the height h of the film. When this theory is implemented for Lennard-Jones liquid films, the new Pi=Pi(h,h(x)) is bounded near a contact line as h-->0. Thus, the singularity in Pi(h) is artificial because it results from extending a theory beyond its range of validity. We also show that the new Pi can capture all three regimes of drop behavior (complete wetting, partial wetting, and pseudo-partial wetting) without altering the signs of the long and short-range interactions. We find that a drop with a precursor film is linearly stable.     

Parameter space minimization methods: applications to Lennard-Jones-dipole-dipole clusters

The morphology of the uniform Lennard-Jones-dipole-dipole cluster with 13 centers (LJDD)13 is investigated over a relatively wide range of values of the dipole moment. We introduce and compare several necessary modifications of the basin-hopping algorithm for global optimization to improve its efficiency. We develop a general algorithm for T=0 Brownian dynamics in curved spaces, and a graph theoretical approach necessary for the elimination of dissociated states. We find that the (LJDD)13 cluster has icosahedral symmetry for small to moderate values of the dipole moment. As the dipole moment increases, however, its morphology shifts to an hexagonal antiprism, and eventually to a ring.     

Combining flavin photocatalysis with parallel synthesis: a general platform to optimize peptides with non-proteinogenic amino acids

Most peptide drugs contain non-proteinogenic amino acids (NPAAs), born out through extensive structure–activity relationship (SAR) studies using solid-phase peptide synthesis (SPPS). Synthetically laborious and expensive to manufacture, NPAAs also can have poor coupling efficiencies allowing only a small fraction to be sampled by conventional SPPS. To gain general access to NPAA-containing peptides, we developed a first-generation platform that merges contemporary flavin photocatalysis with parallel synthesis to simultaneously make, purify, quantify, and even test up to 96 single-NPAA peptide variants via the unique combination of boronic acids and a dehydroalanine residue in a peptide. We showcase the power of our newly minted platform to introduce NPAAs of diverse chemotypes-aliphatic, aromatic, heteroaromatic-directly into peptides, including 15 entirely new residues, and to evolve a simple proteinogenic peptide into an unnatural inhibitor of thrombin by non-classical peptide SAR.

We report a non-classical approach to interrogate peptides with non-proteinogenic amino acids via flavin photocatalysis. We establish a new platform to make, purify, quantify, and biochemically test up to 96 peptide variants in batch.     

Novel method for geometry optimization of molecular clusters: application to benzene clusters

A heuristic and unbiased method for searching optimal geometries of clusters of nonspherical molecules was constructed from the algorithm recently proposed for Lennard-Jones atomic clusters. In the method, global minima are searched by using three operators, interior, surface, and orientation operators. The first operator gives a perturbation on a cluster configuration by moving molecules near the center of mass of a cluster, and the second one modifies a cluster configuration by moving molecules to the most stable positions on the surface of a cluster. The moved molecules are selected by employing a contribution of the molecules to the potential energy of a cluster. The third operator randomly changes the orientations of all molecules. The proposed method was applied to benzene clusters. It was possible to find new global minima for (C6H6)11, (C6H6)14, and (C6H6)15. Global minima for (C6H6)16 to (C6H6)30 are first reported in this article.     

A dynamic lattice searching method with interior operation for unbiased optimization of large Lennard-Jones clusters

For improving the efficiency of dynamic lattice searching (DLS) method for unbiased optimization of large Lennard-Jones (LJ) clusters, a variant of the interior operation (IO) proposed by Takeuchi was combined with DLS. The method is named as DLS-IO. In the method, the IO moves outer atoms with higher energy toward the coordinates center, i.e., (0, 0, 0), of a cluster and a local minimization (LM) follows each IO. This makes the interior atoms more compact and the outer atoms more uniformly distributed with lower potential energy. Therefore, the starting structure for DLS operations is closer to the global optimum compared with the randomly generated structures. On the other hand, a method to identify the central atom is proposed for the central vacancy problem. Optimizations of LJ(500), LJ(561), LJ(660), LJ(665), and LJ(670) were investigated with the DLS-IO, and the structural transition during the optimization was analyzed. It was found that the method is efficient and unbiased for optimization of large LJ clusters, and it may be a promising approach to be universally used for structural optimizations.     

Large-scale parallel configuration interaction. I. Nonrelativistic and scalar-relativistic general active space implementation with application to (Rb-Ba)+

We present a parallel implementation of a string-driven general active space configuration interaction program for nonrelativistic and scalar-relativistic electronic-structure calculations. The code has been modularly incorporated in the DIRAC quantum chemistry program package. The implementation is based on the message passing interface and a distributed data model in order to efficiently exploit key features of various modern computer architectures. We exemplify the nearly linear scalability of our parallel code in large-scale multireference configuration interaction (MRCI) calculations, and we discuss the parallel speedup with respect to machine-dependent aspects. The largest sample MRCI calculation includes 1.5x10(9) Slater determinants. Using the new code we determine for the first time the full short-range electronic potentials and spectroscopic constants for the ground state and for eight low-lying excited states of the weakly bound molecular system (Rb-Ba)+ with the spin-orbit-free Dirac formalism and using extensive uncontracted basis sets. The time required to compute to full convergence these electronic states for (Rb-Ba)+ in a single-point MRCI calculation correlating 18 electrons and using 16 cores was reduced from more than 10 days to less than 1 day.     

Monte Carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical Lennard-Jones fluid

The paper reports a Monte Carlo technique for estimation of the free energies of fluids by sampling on distributions designed for this purpose, rather than on the usual Boltzmann distribution. As an illustration of its use, the free energy of a Lennard-Jones fluid in the liquid-vapour coexistence region has been estimated by relating it to that of the inverse-twelve (soft sphere) fluid, which itself shows no condensation.     

Statistical mechanical theory for the structure of steady state systems: application to a Lennard-Jones fluid with applied temperature gradient

The constrained entropy and probability distribution are given for the structure that develops in response to an applied thermodynamic gradient, as occurs in driven steady state systems. The theory is linear but is applicable to gradients with arbitrary spatial variation. The phase space probability distribution is also given, and it is surprisingly simple with a straightforward physical interpretation. With it, all of the known methods of equilibrium statistical mechanics for inhomogeneous systems may now be applied to determining the structure of nonequilibrium steady state systems. The theory is illustrated by performing Monte Carlo simulations on a Lennard-Jones fluid with externally imposed temperature and chemical potential gradients. The induced energy and density moments are obtained, as well as the moment susceptibilities that give the rate of change of these with imposed gradient and which also give the fluctuations in the moments. It is shown that these moment susceptibilities can be written in terms of bulk susceptibilities and also that the Soret coefficient can be expressed in terms of them.     

First-principle molecular dynamics with ultrasoft pseudopotentials: parallel implementation and application to extended bioinorganic systems

We present a plane-wave ultrasoft pseudopotential implementation of first-principle molecular dynamics, which is well suited to model large molecular systems containing transition metal centers. We describe an efficient strategy for parallelization that includes special features to deal with the augmented charge in the contest of Vanderbilt's ultrasoft pseudopotentials. We also discuss a simple approach to model molecular systems with a net charge and/or large dipole/quadrupole moments. We present test applications to manganese and iron porphyrins representative of a large class of biologically relevant metalorganic systems. Our results show that accurate density-functional theory calculations on systems with several hundred atoms are feasible with access to moderate computational resources.     

Larger water clusters with edges and corners on their way to ice: structural trends elucidated with an improved parallel evolutionary algorithm

For the difficult task of finding global minimum energy structures for molecular clusters of nontrivial size, we present a highly efficient parallel implementation of an evolutionary algorithm. By completely abandoning the traditional concept of generations and by replacing it with a less rigid pool concept, we have managed to eliminate serial bottlenecks completely and can operate the algorithm efficiently on an arbitrary number of parallel processes. Nevertheless, our new algorithm still realizes all of the main features of our old, successful implementation. First tests of the new algorithm are shown for the highly demanding problem of water clusters modeled by a potential with flexible, polarizable monomers (TTM2-F). For this problem, our new algorithm not only reproduces all of the global minima proposed previously in considerably less CPU time but also leads to improved proposals in several cases. These, in turn, qualitatively change our earlier predictions concerning the transitions from all-surface structures to cages with a single interior molecule, and from one to two interior molecules. Furthermore, we compare preliminary results up to n = 105 with locally optimized cuts from several ice modifications. This comparison indicates that relaxed ice structures may start to be competitive already at cluster sizes above n = 90.     

Many-body interaction analysis: algorithm development and application to large molecular clusters

A completely automated algorithm for performing many-body interaction energy analysis of clusters (MBAC) [M. J. Elrodt and R. J. Saykally, Chem. Rev. 94, 1975 (1994); S. S. Xantheas, J. Chem. Phys. 104, 8821 (1996)] at restricted Hartree-Fock (RHF)/MA Plesset 2nd order perturbation theory (MP2)/density functional theory (DFT) level of theory is reported. Use of superior guess density matrices (DM's) for smaller fragments generated from DM of the parent system and elimination of energetically insignificant higher-body combinations, leads to a more efficient performance (speed-up up to 2) compared to the conventional procedure. MBAC approach has been tested out on several large-sized weakly bound molecular clusters such as (H(2)O)(n), n=8, 12, 16, 20 and hydrated clusters of amides and aldehydes. The MBAC results indicate that the amides interact more strongly with water than aldehydes in these clusters. It also reconfirms minimization of the basis set superposition error for large cluster on using superior quality basis set. In case of larger weakly bound clusters, the contributions higher than four body are found to be repulsive in nature and smaller in magnitude. The reason for this may be attributed to the increased random orientations of the interacting molecules separated from each other by large distances.     

Probing isomer interconversion in anionic water clusters using an Ar-mediated pump-probe approach: Combining vibrational predissociation and velocity-map photoelectron imaging spectroscopies

We present the first results from an experiment designed to explore barriers for interconversion between isomers of cluster anions using an Ar-cluster mediated pump-probe technique. In this approach, anions are generated with many Ar atoms attached, and one of the isomers present is selectively excited by tuning an infrared laser to one of the isomer's characteristic vibrational resonances. The excited cluster is then cooled by evaporation of Ar atoms, and the isomer distribution in the lighter daughter ions is measured after secondary mass selection by recording their photoelectron spectra using velocity-map imaging. We apply the method to the water hexamer anion, (H(2)O)(6) (-), which is known to occur in two isomeric forms with different electron-binding energies. We find that conversion of the high-binding (type I) form to the low-binding (type II) isomer is not efficiently driven in (H(2)O)(6) (-) with excitation energies in the 0.4 eV range even though it is possible to create both isomers in abundance in the ion source. This observation is discussed in the context of the competition between isomerization and electron autodetachment, which depends on the relative positions of the neutral and ionic potential surfaces along the isomerization pathway. Application of the method to the more complex heptamer ion, however, does reveal that interconversion is available among the highest binding isomer classes (I and I(')).     

Theory of electrode polarization: application to parallel plate cell dielectric spectroscopy experiments

An extension of the model for electrode polarization of Cirkel et al. [Physica A 235 (1997) 269] is given. The problem is solved using both classical boundary conditions and the new boundary conditions using excess densities presented in a previous paper [J. Phys. Chem. B 105 (2001) 11743]. In the present paper, the electrodes are supposed to be ideal, meaning that charge transfer or adsorption are not considered. The advantage of the new boundary conditions lies in the possibility to extend to more complicated situations including for instance specific ion adsorption. We prove that the new boundary conditions and classical ones give the same results. A comparison of the model predictions, involving no adjustable parameters, experimental dielectric spectroscopy data is performed and fairly good agreement is found.     

Communication: The application of the global isomorphism to the study of liquid-vapor equilibrium in two and three-dimensional Lennard-Jones fluids

We analyze the interrelation between the coexistence curve of the Lennard-Jones fluid and the Ising model in two and three dimensions within the global isomorphism approach proposed earlier [V. L. Kulinskii, J. Phys. Chem. B 114, 2852 (2010)]. In case of two dimensions, we use the exact Onsager result to construct the binodal of the corresponding Lennard-Jones fluid and compare it with the results of the simulations. In the three-dimensional case, we use available numerical results for the Ising model for the corresponding mapping. The possibility to observe the singularity of the binodal diameter is discussed.