Practical quantum mechanics-based fragment methods for predicting molecular crystal properties

Significant advances in fragment-based electronic structure methods have created a real alternative to force-field and density functional techniques in condensed-phase problems such as molecular crystals. This perspective article highlights some of the important challenges in modeling molecular crystals and discusses techniques for addressing them. First, we survey recent developments in fragment-based methods for molecular crystals. Second, we use examples from our own recent research on a fragment-based QM/MM method, the hybrid many-body interaction (HMBI) model, to analyze the physical requirements for a practical and effective molecular crystal model chemistry. We demonstrate that it is possible to predict molecular crystal lattice energies to within a couple kJ mol(-1) and lattice parameters to within a few percent in small-molecule crystals. Fragment methods provide a systematically improvable approach to making predictions in the condensed phase, which is critical to making robust predictions regarding the subtle energy differences found in molecular crystals.     

Anisotropic interfacial free energies of the hard-sphere crystal-melt interfaces

We present a reliable method to define the interfacial particles for determining the crystal-melt interface position, which is the key step for the crystal-melt interfacial free energy calculations using capillary wave approach. Using this method, we have calculated the free energies gamma of the fcc crystal-melt interfaces for the hard-sphere system as a function of crystal orientations by examining the height fluctuations of the interface using Monte Carlo simulations. We find that the average interfacial free energy gamma(0) = 0.62 +/- 0.02k(B)T/sigma(2) and the anisotropy of the interfacial free energies are weak, gamma(100) = 0.64 +/- 0.02, gamma(110) = 0.62 +/- 0.02, gamma(111) = 0.61 +/- 0.02k(B)T/sigma(2). The results are in good agreement with previous simulation results based on the calculations of the reversible work required to create the interfaces (Davidchack and Laird, Phys. Rev. Lett. 2000, 85, 4571). In addition, our results indicate gamma(100) > gamma(110) > gamma(111) for the hard-sphere system, similar to the results of the Lennard-Jones system.     

Thermodynamic analysis of pure and impurity doped pentaerythritol tetranitrate crystals grown at room temperature

Pentaerythritol tetranitrate (PETN) powders are used to initiate other explosives. During long-term storage, changes in powder properties can cause changes in the initiation performance. Changes in the morphology and surface area of aging powders are observed due to sublimation and growth of PETN crystals through coarsening mechanisms, (e.g. Ostwald ripening, sintering, etc.). In order to alleviate the sublimation of PETN crystals under service conditions, stabilization methods such as thermal cycling and doping with certain impurities during or after the crystallization of PETN have been proposed. In this report we present our work on the effect of impurities on the morphology and activation energy of the PETN crystals. The pure and impurity doped crystals of PETN were grown from supersaturated acetone solution by solvent evaporation technique at room temperature. The difference in the morphology of the impurity-doped PETN crystal compared to pure crystal was examined by optical microscopy. The changes in the activation energies and the evaporation rates are determined by thermogravimetry (TG). Our activation energies of evaporation agree with earlier reported enthalpies of vaporization. The morphology and activation energy of PETN crystals doped with Ca, Na, and Fe cations are similar to that for pure PETN crystal, whereas the Zn-ion-doped PETN crystals have different morphology and decreased activation energy.     

The Role of Liquid Water in Crystalline Hydrate Dehydration: Copper sulphate pentahydrate

Simultaneous Differential Thermal Analysis/Thermogravimetric experiments carried out on one large single crystal, several small single crystals and powdered crystals of pentahydrate copper sulphate have been used to demonstrate the role that retained liquid water plays in maintaining crystal morphology during dehydration. Measured activation energies for stepwise dehydration in the system show the presence of solution-based transformations provide lower energy paths for the dehydration steps and stress relieving mechanisms. Skeletal anhydrous crystals from large-sized pentahydrate copper sulphate have the same morphology as the starting crystal on complete dehydration at controlled heating rates as long as a solution phase is maintained within the crystal during decomposition. The athermal activation energies, in this work, are in agreement with those obtained by isothermal methods as long as coincident reaction paths for the two techniques are maintained. The literature has been reviewed in the light of this work and a three-stage process is presented to rationalise the conflicting information obtained by workers using a variety of different experimental techniques.This revised version was published online in November 2005 with corrections to the Cover Date.     

New simulation model of multicomponent crystal growth and inhibition

We review a novel computational model for the study of crystal structures both on their own and in conjunction with inhibitor molecules. The model advances existing Monte Carlo (MC) simulation techniques by extending them from modeling 3D crystal surface patches to modeling entire 3D crystals, and by including the use of "complex" multicomponent molecules within the simulations. These advances makes it possible to incorporate the 3D shape and non-uniform surface properties of inhibitors into simulations, and to study what effect these inhibitor properties have on the growth of whole crystals containing up to tens of millions of molecules. The application of this extended MC model to the study of antifreeze proteins (AFPs) and their effects on ice formation is reported, including the success of the technique in achieving AFP-induced ice-growth inhibition with concurrent changes to ice morphology that mimic experimental results. Simulations of ice-growth inhibition suggest that the degree of inhibition afforded by an AFP is a function of its ice-binding position relative to the underlying anisotropic growth pattern of ice. This extended MC technique is applicable to other crystal and crystal-inhibitor systems, including more complex crystal systems such as clathrates.     

A new model for simulating 3-d crystal growth and its application to the study of antifreeze proteins

A novel computational technique for modeling crystal formation has been developed that combines three-dimensional (3-D) molecular representation and detailed energetics calculations of molecular mechanics techniques with the less-sophisticated probabilistic approach used by statistical techniques to study systems containing millions of molecules undergoing billions of interactions. Because our model incorporates both the structure of and the interaction energies between participating molecules, it enables the 3-D shape and surface properties of these molecules to directly affect crystal formation. This increase in model complexity has been achieved while simultaneously increasing the number of molecules in simulations by several orders of magnitude over previous statistical models. We have applied this technique to study the inhibitory effects of antifreeze proteins (AFPs) on ice-crystal formation. Modeling involving both fish and insect AFPs has produced results consistent with experimental observations, including the replication of ice-etching patterns, ice-growth inhibition, and specific AFP-induced ice morphologies. Our work suggests that the degree of AFP activity results more from AFP ice-binding orientation than from AFP ice-binding strength. This technique could readily be adapted to study other crystal and crystal inhibitor systems, or to study other noncrystal systems that exhibit regularity in the structuring of their component molecules, such as those associated with the new nanotechnologies.     

"Exact" surface free energies of iron surfaces using a modified embedded atom method potential and lambda integration

Previously a new universal lambda-integration path and associated methodology was developed for the calculation of "exact" surface and interfacial free energies of solids. Such a method is in principle applicable to any intermolecular potential function, including those based on ab initio methods, but in previous work the method was only tested using a relatively simple embedded atom method iron potential. In this present work we apply the new methodology to the more sophisticated and more accurate modified embedded atom method (MEAM) iron potential, where application of other free- energy methods would be extremely difficult due to the complex many-body nature of the potential. We demonstrate that the new technique simplifies the process of obtaining "exact" surface free energies by calculating the complete set of these properties for the low index surface faces of bcc and fcc solid iron structures. By combining these data with further calculations of liquid surface tensions we obtain the first complete set of exact surface free energies for the solid and liquid phases of a realistic MEAM model system. We compare these predictions to various experimental and theoretical results.     

表面能与晶体生长/溶解动力学研究的新动向

界面现象使物质在结晶过程中出现了临界现象.但最近的研究指出在物质溶解过程中,在表面能量的控制下也存在着临界现象以及尺寸效应.实验发现,当晶体自身小到一 定的程度时(通常在纳米尺度上并和临界蚀坑的大小相近),在溶解过程中其速度会自发降 低,反应被抑制乃至停止.尽管在热力学上表面能的因素可以赋予小颗粒晶体较大的溶解度 ,但表面能却也能通过对临界条件的控制而使这些微粒在动力学上不被溶解.这个发现不仅 解决了纳米颗粒在水溶液中稳定性的问题,而且还从动力学的角度解释了生物矿物选择纳米 尺度作为其基本构成单元的原因.由于表面能和晶体生长/溶解的动力学有着密切的关系, 我们可以通过对表面能的调节来修改它们的动力学速度和晶体的形貌.反过来,也可以用动力学的方法来测定表面能及表面吸附/脱附常数等.相对于常规的界面研究手段,通过生 长和溶解动力学途径所得的数据有着很好的可靠性及重复性.我们认为,晶体生长和溶解的 动力学和表面能的研究相结合,不仅为界面研究提供了新的思路和方法,而且也会推动晶体生长和材料科学的发展.     

Properties of inhibitors of methane hydrate formation via molecular dynamics simulations

Within the framework of a proposed two-step mechanism for hydrate inhibition, the energy of binding of four inhibitor molecules (PEO, PVP, PVCap, and VIMA) to a hydrate surface is estimated with molecular dynamic simulations. One key feature of this proposed mechanism is that the binding of an inhibitor molecule to the surface of an ensuing hydrate crystal disrupts growth and therein crystallization. It is found through the molecular dynamic simulations that inhibitor molecules that experimentally exhibit better inhibition strength also have higher free energies of binding, an indirect confirmation of our proposed mechanism. Inhibitors increasing in effectiveness, PEO < PVP < PVCap < VIMA, have increasingly negative (exothermic) binding energies of -0.2 < -20.6 < -37.5 < -45.8 kcal/mol and binding free energies of increasing favorability (+0.4 approximately = +0.5 < -9.4 < -15.1 kcal/mol). Furthermore, the effect of an inhibitor molecule on the local liquid water structure under hydrate-forming conditions was examined and correlated to the experimental effectiveness of the inhibitors. Two molecular characteristics that lead to strongly binding inhibitors were found: (1) a charge distribution on the edge of the inhibitor that mimics the charge separation in the water molecules on the surface of the hydrate and (2) the congruence of the size of the inhibitor with respect to the available space at the hydrate-surface binding site. Equipped with this molecular-level understanding of the process of hydrate inhibition via low-dosage kinetic hydrate inhibitors we can design new, more effective inhibitor molecules.     

Phase field theory of interfaces and crystal nucleation in a eutectic system of fcc structure: I. Transitions in the one-phase liquid region

The phase field theory (PFT) has been applied to predict equilibrium interfacial properties and nucleation barrier in the binary eutectic system Ag-Cu using double well and interpolation functions deduced from a Ginzburg-Landau expansion that considers fcc (face centered cubic) crystal symmetries. The temperature and composition dependent free energies of the liquid and solid phases are taken from CALculation of PHAse Diagrams-type calculations. The model parameters of PFT are fixed so as to recover an interface thickness of approximately 1 nm from molecular dynamics simulations and the interfacial free energies from the experimental dihedral angles available for the pure components. A nontrivial temperature and composition dependence for the equilibrium interfacial free energy is observed. Mapping the possible nucleation pathways, we find that the Ag and Cu rich critical fluctuations compete against each other in the neighborhood of the eutectic composition. The Tolman length is positive and shows a maximum as a function of undercooling. The PFT predictions for the critical undercooling are found to be consistent with experimental results. These results support the view that heterogeneous nucleation took place in the undercooling experiments available at present. We also present calculations using the classical droplet model [classical nucleation theory (CNT)] and a phenomenological diffuse interface theory (DIT). While the predictions of the CNT with a purely entropic interfacial free energy underestimate the critical undercooling, the DIT results appear to be in a reasonable agreement with the PFT predictions.     

High-resolution electron energy loss spectroscopy study of O-Cu(410)

We have investigated oxygen adsorption on Cu(410) by high-resolution electron energy loss spectroscopy, dosing O2 with a supersonic molecular beam at different surface temperatures and for different angles of incidence and beam energies or by backfilling. In the investigated crystal temperature range (127 < T < 570 K), adsorption is always dissociative. Depending on T, impact energy, and angle of incidence, the oxygen atoms end up in different adsorption configurations, characterized by different vibrational signatures. In particular, at grazing incidence when only the step edge is exposed to O2, the adatoms end up initially preferentially at the step edge. An ordered overlayer forms at half monolayer coverage when the adsorbate is mobile. Oxide patches develop eventually for large exposures performed by backfilling and at high crystal temperature.     

Calculations of crystal-melt interfacial free energies by nonequilibrium work measurements

We developed a multistep thermodynamic perturbation method to compute the interfacial free energies by nonequilibrium work measurements with cleaving potential procedure. Using this method, we calculated the interfacial free energies of different crystal orientations for the Lennard-Jones system. Our results are in good agreement with the results by thermodynamic integration method. Compared with thermodynamic integration method, the multistep thermodynamic perturbation method is more efficient. For each stage of the cleaving process, only a few thermodynamic perturbation steps are needed, and there is no requirement on the reversibility of the path.     

Effect of bicarbonate ions on the crystallization of calcite on self-assembled monolayers

We use molecular dynamics simulations to investigate the nucleation of calcite crystals on self-assembled monolayers. We show how the presence of bicarbonate ions adsorbed on the monolayer surface can both aid nucleation and control the orientation of the growth of the crystal. Using a simple model of the nucleation process and calculated interfacial energies, we calculate the enhancement (with respect to the homogeneous nucleation rate) of the nucleation of calcite on the (012) and (0001) faces. The calculations show clearly that the (012) face is favored over the (0001) face and that the nucleation rate is enhanced for self-assembled monolayers made from molecules containing an even number of carbon atoms in the alkyl chain over those containing an odd number.     

Thermal and spectral properties and induction period, interfacial energy and nucleation parameters of solution grown anthracene

Anthracene crystals were grown by solution growth technique by adopting slow evaporation method from the solvents CS₂, CCl₄ and CHCl₃. The induction period was measured at various super saturations, and hence the interfacial energies were evaluated. Using the interfacial tension value, the nucleation parameters such as radius of the critical nuclei (r*), the Gibbs free energy change for the formation of a critical nucleus (?G*) and the number of molecules in the critical nucleus (i*) were also calculated for all these solvents at two different temperatures. The effect of surface tension, viscosity and density of these solvents are correlated with interfacial tension. The solution grown crystals were subjected to UV, FTIR, NMR and X-ray diffraction studies. The purity and high-thermal stability of the grown crystals were determined using thermal analysis.     

Calculation of solid-liquid interfacial free energy: a classical nucleation theory based approach

We present a simple approach to calculate the solid-liquid interfacial free energy. This new method is based on the classical nucleation theory. Using the molecular dynamics simulation, we employ spherical crystal nuclei embedded in the supercooled liquids to create an ideal model of a homogeneous nucleation. The interfacial free energy is extracted by fitting the relation between the critical nucleus size and the reciprocal of the critical undercooling temperature. The orientationally averaged interfacial free energy is found to be 0.302+/-0.002 (in standard LJ unit). The temperature dependence of the interfacial free energy is also obtained in this work. We find that the interfacial free energy increases slightly with increasing temperature. The positive temperature coefficient of the interfacial free energy is in qualitative agreement with Spaepen's analysis [Solid State Phys. 47, FS181 (1994)] and Turnbull's empirical estimation [J. Appl. Phys. 21, 1022 (1950)].     

XDM-corrected hybrid DFT with numerical atomic orbitals predicts molecular crystal lattice energies with unprecedented accuracy

Molecular crystals are important for many applications, including energetic materials, organic semiconductors, and the development and commercialization of pharmaceuticals. The exchange-hole dipole moment (XDM) dispersion model has shown good performance in the calculation of relative and absolute lattice energies of molecular crystals, although it has traditionally been applied in combination with plane-wave/pseudopotential approaches. This has limited XDM to use with semilocal functional approximations, which suffer from delocalization error and poor quality conformational energies, and to systems with a few hundreds of atoms at most due to unfavorable scaling. In this work, we combine XDM with numerical atomic orbitals, which enable the efficient use of XDM-corrected hybrid functionals for molecular crystals. We test the new XDM-corrected functionals for their ability to predict the lattice energies of molecular crystals for the X23 set and 13 ice phases, the latter being a particularly stringent test. A composite approach using a XDM-corrected, 25% hybrid functional based on B86bPBE achieves a mean absolute error of 0.48 kcal mol⁻¹ per molecule for the X23 set and 0.19 kcal mol⁻¹ for the total lattice energies of the ice phases, compared to recent diffusion Monte-Carlo data. These results make the new XDM-corrected hybrids not only far more computationally efficient than previous XDM implementations, but also the most accurate density-functional methods for molecular crystal lattice energies to date.

A new implementation of the XDM dispersion model within FHI-aims allows efficient use of XDM-corrected hybrid functionals, found to be the most accurate DFT methods for molecular crystal lattice energies to date.     

Mechanical model for anisotropic curved interfaces with applications to surfactant-laden liquid-liquid crystal interfaces

A mechanical model for anisotropic curved interfaces, applicable to thermodynamically closed surfactant-laden liquid-liquid crystal interfaces is developed. The model takes into account the mechanical effects due to surface bending and surface tilting (anchoring) and incorporates liquid crystal anisotropy into classical fluid membrane mechanics. In the absence of the aligned liquid crystal, the model converges to the fluid membrane mechanical model, and in the absence of surfactant, it converges to the nematic interface mechanical model. Use of the well-known Helfrich-Rapini-Papoular surface energies leads to the Laplace equation for anisotropic curved interfaces, whose material limits are the vesicle shape equation and the liquid crystal Herring equation. Applications of the model to shape selection in liquid drops embedded in aligned nematic liquid crystals illustrates the competition between anchoring and bending and shows how anisotropic surface tension distorts the droplet and how bending tends to restore the spherical shape. This theory presented in this article shows that the interaction of interfacial anchoring and bending creates new regimes in classical fluid membrane mechanics.     

Direct calculation of the crystal-melt interfacial free energy via molecular dynamics computer simulation

We review our recent work on the direct calculation of the interfacial free energy, gamma, of the crystal-melt interface via molecular dynamics computer simulation for a number of model systems. The value of gamma as a function of crystal orientation is determined using a thermodynamic integration technique employing moving cleaving walls [Phys. Rev. Lett. 2000, 85, 4751]. The calculation is sufficiently precise to resolve the small anisotropy in gamma, which is crucial in determining the kinetics and morphology of dendritic growth. We report values of gamma for the hard-sphere and Lennard-Jones systems, as well as recent results on the series of inverse-power potentials. For the inverse sixth-, seventh-, and eighth-power systems, we determine gamma for both fcc and bcc crystal structures. For these systems, the bcc-melt gamma is lower than that for fcc crystals by about 25%, consistent with recent experiments and computer simulations on fcc-forming systems that show preferential formation of bcc nuclei in the initial stages of crystallization. In addition, we show that our results give a molecular interpretation of Turnbull's rule, which is the empirical relationship between gamma and the enthalpy of fusion.     

Interfacial properties of statistical copolymer brushes in contact with homopolymer melts

We use polymer self-consistent field theory to quantify the interfacial properties of random copolymer brushes (AB) in contact with a homopolymer melt chemically identical to one of the blocks (A). We calculate the interfacial widths and interfacial energies between the melt and the brush as a function of the relative chain sizes, grafting densities, compositions of the random copolymer in the brush, and degree of chemical incompatibility between the A and B species. Our results indicate that the interfacial energies between the melt and the brush increase (signifying expulsion of the free chains from the brush) with increasing grafting density, chemical incompatibility between A and B components, and size of the free chains relative to the grafted chains. We also compare the interfacial energies of random copolymers of different sequence characteristics and find that, except for the case of very blocky or proteinlike chains, blockiness of the copolymer has only little effect on interfacial properties. Our results for interfacial energies are rationalized based on the concept of an "effective volume fraction" of the brush copolymers, f(eff), which quantifies the chemical composition of the brush segments in the interfacial zone between the brush and melt copolymers. Using this concept, we modify the strong-stretching theory of brush-melt interfaces to arrive at a simple model whose results qualitatively agree with our results from self-consistent field theory. We discuss the ramifications of our results for the design of neutral surfaces.