Diabat Interpolation for Polymorph Free-Energy Differences

Subscriber access provided by Iowa State University | Library

Letter

Diabat Interpolation for Polymorph Free Energy Differences Kartik Kamat, and Baron Peters J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.6b02795 • Publication Date (Web): 17 Jan 2017 Downloaded from http://pubs.acs.org on January 23, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Diabat Interpolation for Polymorph Free Energy Differences Kartik Kamat† and Baron Peters∗,†,‡ †Department of Chemical Engineering, University of California Santa Barbara, CA 93106 ‡Department of Chemistry and Biochemistry, University of California Santa Barbara, CA 93106 E-mail: [email protected]

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Existing methods to compute free energy differences between polymorphs use harmonic approximations, advanced non-Boltzmann bias sampling techniques, and/or multistage free energy perturbations. This work demonstrates how Bennett’s diabat interpolation method (J. Comp. Phys. 1976, 22, 245) can be combined with energy gaps from lattice switch Monte Carlo techniques (Phys. Rev. E 2001, 61, 906) to swiftly estimate polymorph free energy differences. The new method requires only two unbiased molecular dynamics simulations, one for each polymorph. To illustrate the new method, we compute the free energy difference between face centered cubic (FCC) and body centered cubic (BCC) polymorphs for a Gaussian core solid. We discuss the justification for parabolic models of the free energy diabats and similarities to methods that have been used in studies of electron transfer.

Graphical TOC Entry


Page 2 of 20

Page 3 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Polymorphs are different crystal structures for the same molecule or compound. 1 Different polymorphs can exhibit different mechanical properties, 2,3 different optical properties, 4,5 different solubilities, 6,7 different melting temperatures, 8,9 different catalytic activities, 10,11 different crystal growth habits, 12,13 etc. Because it affects so many properties, polymorphism is important to pharmaceuticals, 14–17 biomineralization, 18,19 metallurgy, 20,21 geology, 22,23 etc. The key challenges in computational studies of polymorphism are: 1) discovery/prediction of new polymorphs, 24–28 2) quantifying the free energy differences between polymorphs, 29–33 and 3) predicting nucleation and growth kinetics for each polymorph. 34,35 Significant progress has been made in polymorph prediction and some progress has been made in crystallization kinetics, but successful kinetic predictions have been limited to model systems and/or to simple molecules with accurate forcefields and precisely quantified driving forces. Typical free energy differences between observed polymorphs are only a couple of kJ/mol, 36 so subtle factors like dispersion, 37 polarization, 31 and anharmonic vibrational effects 38 are important. For more complex molecules, the difficulty of computing free energy differences between polymorphs hampers computational studies of nucleation and growth kinetics. Reviews on methods to compute polymorph free energy differences can be found elsewhere, 39–41 but the most prevalent techniques are harmonic approximations, 42 quasi-harmonic expansions, 43–45 and fully anharmonic calculations. The anharmonic and quasi-anharmonic calculations can be performed with accurate ab initio Hamiltonians. 46 This paper focuses on methods that can treat fully anharmonic modes, but which usually require force fields to accelerate sampling. Among these, the Frenkel-Ladd approach 29 and the related Einstein molecule approach 47 compute the reversible work to transform polymorph a to an Einstein crystal and then to polymorph b. They obtain the small free energy difference between polymorphs from the difference between two larger free energies along long paths to and from the Einstein crystal reference system. Tan et al. 32 developed harmonically targeted temperature perturbation (HTTP) methods to compute (classical) absolute free energies of polymorphs from a series of perturbations starting at very low temperature where the har-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 20

monic approximation is valid. Martonak et al. 30 and Tuckerman et al. 48 have used lattice vectors as order parameters with biased sampling techniques to sample a reversible (but dynamically fictitious) path from one polymorph to another. Each of the above methods can provide highly accurate free energy differences between polymorphs, but they are rather demanding calculations. Here we note that lattice switch moves adapted from the Lattice-Switch Monte Carlo (LSMC) method of Bruce et al. 33,49 can be used in a very different, and comparatively easy, route to the free energy difference between polymorphs. First, we briefly describe the LSMC moves. The complete position vector of a ~ iα ) and a displacement vector particle (~riα ) can be written as a sum of a reference position (R (~ui ). ~ α + ~ui ~riα = R i

(1)

Here α refers to a particular phase. A lattice switch move attempts to swap the reference ~ α without changing the displacements. The acceptance probability of an lattice positions R i LSMC move depends on the energy gap

∆E(~ u) = Ea (~ u) − Eb (~ u)

(2)

i.e the potential energy difference between polymorphs a and b with the same displacement vector u ~ in each polymorph. The typical displacements of atoms during an unbiased simulation of polymorph a will be different from the typical displacements in polymorph b, so that LSMC moves have large energy gaps and extremely low acceptance probabilities. Jackson et al. demonstrated that LSMC simulations can be biased along ∆E to generate the “adiabatic” Landau free energy FL (∆E). 50 Biased LSMC free energy calculations have been performed for core-softened fluids, 51 Lennard-Jones solids, 50 solid phases of butane 52 and BCC and HCP phases of zirconium. 53 Like other methods for computing free energy differences between polymorphs, the LSMC path is a non-physical route. However it suggests interesting analogies to methods that use


Page 5 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


energy gaps to compute free energy difference in electron transfer applications. 54–56 In an electron transfer reaction, 57 the electron jumps from the donor atom to an acceptor atom with all atoms at fixed positions. The LSMC move analogously jumps between the reference lattices while maintaining fixed displacements (~ u). Motivated by this analogy, let ρa (∆E) and ρb (∆E) be distributions of the energy gap that are generated during unbiased simulations in states a and b, respectively. The well-known Bennett Acceptance Ratio (BAR) method 58 requires significant overlap between ρa (∆E) and ρb (∆E), and these will rarely overlap with each other. In the same paper, Bennett suggests an alternative diabat interpolation scheme for cases where ρa (∆E) and ρb (∆E) are non-overlapping on the ∆E axis, but sufficiently smooth to permit interpolation. 58 A theoretically motivated model for the diabats helps to bridge the gap between unbiased ρa (∆E) and ρb (∆E) histograms. The principal advantages of diabat interpolation are its low computational cost and its easy implementation. It requires just two unconstrained simulations (one for each state), while other methods require multistage perturbations or advanced sampling techniques. Diabat interpolation has been used in electron transfer studies, e.g. Blumberger et al. 54 used parabolic diabat models to interpolate the Marcus curves for electron transfer between Ag + /Ag 2+ ions in water. In Bennett’s formulation of the diabat interpolation approach, ρa (∆E) and ρb (∆E) are distributions of the energy gap from two unbiased simulations with different Hamiltonians. In the proposed extension to polymorphism, ρa (∆E) and ρb (∆E) are distributions of the energy gap from unbiased simulations with the same Hamiltonian, but with two different reference lattices. Despite this difference the diabat interpolation framework is applicable because: (i) the basins corresponding to states a and b are distant and non-overlapping parts of configurational space, (ii) the potential energy surfaces in the neighborhood of reference lattices a and b are different, and (iii) there is a one-to-one mapping (with unit Jacobian) between configurations in states a and b. For all practical purposes, the lattice switch can be viewed as a change in the Hamiltonian.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 20

The diabatic free energy curves for two different polymorphs in the NVT ensemble can be written as: Z βFa (∆E) = −ln

d~ u e−βEa (~u) δ(∆E(~ u) − ∆E)

(3)

d~ u e−βEb (~u) δ(∆E(~ u) − ∆E)

(4)

a

Z βFb (∆E) = −ln b

Here β = 1/kB T , a and b correspond to the different polymorphs, and the

R a

and

R b

indicate

~ a and R ~ b . Subtracting the integration over displacements applied to reference lattices R Equation 4 from Equation 3 yields: 58

Fa (∆E) − Fb (∆E) = kB T ln

heβ∆E(~u) δ(∆E(~ u) − ∆E)ia hδ(∆E(~ u) − ∆E)ia

(5)

where h ia indicates an average over the ensemble of equilibrium displacements in state a. Equation 5 can be further simplified to give :

Fa (∆E) − Fb (∆E) = ∆E

(6)

From Equation 6 it can be inferred that the two diabats must intersect each other at ∆E = 0. Although Fb (∆E) and Fa (∆E) have minima at different ∆E values, the vertical difference between them is always given by the energy gap itself and this dramatically facilitates the polymorph free energy calculation. An unbiased simulation in one basin, say phase a, directly generates equilibrium data on the Fa (∆E) curve, but because of Equation 6, the same simulation data provides a nonequilibrium section of the Fb (∆E) curve. Similar piecewise portions of the two diabats can be obtained from an unbiased simulation of phase b. Thus two unbiased simulations yield four piecewise sections of the free energy diabats Fa (∆E) and Fb (∆E). The two non-equilibrium sections would have been very difficult to aquire by other means. A parabolic model for the diabats is motivated by the Central Limit Theorem. 59 The energy gap is the difference between two random variables Ea (~ u) and Eb (~ u) which are sums over thousands of atoms 6 ACS Paragon Plus Environment

Page 7 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that have similar interactions with their neighbors. As such, we expect Ea and Eb to both be approximately Gaussian random variables, and therefore the energy gap should also be a Gaussian random variable. 59 Based on the above considerations we propose parabolic models for the diabatic free energy curves: Fa (∆E) = λ(∆E)2 + φ∆E + F0

(7)

Fb (∆E) = λ(∆E)2 + (φ − 1)∆E + F0

(8)

and

The parabolic model for the diabats facilitates interpolation across the gap between typical values of ∆E in states a and b. When the diabats are not perfect parabolas, the deviation can be modeled by replacing λ in Equations 7 and 8 with λ0 (1 + ε∆E). The non-parabolic corrections will be important when ρa (∆E) and ρb (∆E) have different widths. The difference is not pronounced in this work, but adiabatic LSMC results for the BCC and HCP polymorphs of zirconium 53 provide an example where the non-parabolic corrections would be needed. Figure 1 illustrates the analogy between the diabats in electron transfer and the diabats that we use to estimate polymorph free energy differences. Figure 1 also shows the two electronic configurations and fixed atom positions during electron transfer, as well as the two reference lattices with fixed displacements during a lattice switch. Both procedures convert a typical equilibrium reactant (product) configuration into a non-equilibrium product (reactant) configuration. To illustrate the diabat interpolation approach, we compute the free energy difference between the face centered cubic (FCC) and body centered cubic (BCC) phases for Gaussian core solids. The Gaussian core potential 60 is :

U (r) = e−r

2 /σ 2


(9)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: A) Electron transfer from a typical configuration of the donor D− , acceptor A, and solvent incurs a large energy gap. (B) Parabolic free energy diabats with a typical equilibrium location in the reactant state marked by 5. The corresponding non-equilibrium position on the product diabat is also indicated. The reactant and product state diabats cross each other at ∆E = 0. (C) Lattice switch moves from equilibrium configurations in polymorph a map typical equilibrium displacements to high energy configurations in state b. The circled blue atoms (l) show how typical displacements can lead to a large energy gap. where determines the repulsion strength and σ determines the repulsion range. The numbering of atoms in the FCC and the BCC phases determines the one-to-one mapping between displacements during a lattice switch. Different numbering schemes can potentially influence the distributions of energy gaps and the efficiency of the interpolation procedure. The importance of the numbering scheme for efficiency was also reported in regular LSMC simulations by Bridgwater and Quigley. 52 Each permutation of the atom labels creates a symmetry-equivalent version of the system. Accordingly, each polymorph can be represented by any one of several symmetry equivalent minima and basins. The schematic below depicts the basins corresponding to two polymorphs a and b, and their symmetric representations a0 and b0 that result from switching the atom labels. Note that the stiff/soft modes in each basin are oriented differently from those in its symmetric labelswapped representation. We are free to choose the atom labels so that state a maps to b or to the alternative representation b0 . The atom labels in the schematic are not ideal for state 8 ACS Paragon Plus Environment

Page 8 of 20

Page 9 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Illustrating required symmetries on the energy landscape for polymorphs a, b, and their symmetric representations a0 and b0 with permuted atom labels. The mapping from polymorph a directly to b is not ideal because the distribution of displacements u (red) will have very different shapes in a and b. Instead, we can relabel the atoms in state b so that a lattice switch maps polymorph a directly to the b0 representation of polymorph b. b, because a lattice switch from a to b maps the softest displacement direction in a onto the stiffest displacement direction in b. By switching the atom labels so that a maps to b0 , displacements along the soft modes of a map to displacements along the soft modes of b0 . The remainder of this letter refers only to a and b. Alternative representations need not be considered once an appropriate mapping has been adopted. We note that the optimal atom-to-atom mapping might be identified from the following quantitative considerations. The energy gap can be expanded in the displacement u ~ as 1 † ~ (Ha − Hb )~ u + ... ∆E(~ u) = ∆E(~0) + u 2

(10)

~ a and R ~ b , respectively. Each where Ha and Hb are the force constant matrices at locations R relabeling of the atoms in b corresponds to a reordering of the rows and columns of matrix Hb . The ideal labels should make Hb as close as possible to Ha . Additional work is needed to identify the appropriate matrix distance metric for minimizing the difference between Ha and Hb over all possible mappings. Here we represent the FCC phase with a special distorted



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

version of the BCC unit cell as shown in Figure 3B. The unconventional unit cell creates

Figure 3: A) Two conventional face centered cubic (FCC) unit cells where the face centered atoms on the front and back have been neglected for clarity. B) The figure illustrates an alternative FCC unit cell that resembles a distorted BCC unit cell. C) A body centered cubic (BCC) unit cell. a one-to-one mapping between the FCC and BCC phases such that neighboring atoms in polymorph a map to neighboring atoms in polymorph b. Simulations of the Gaussian core particles were carried out in the NVT (canonical) ensemble. The temperature was controlled with the Nose-Hoover 61,62 style thermostat as implemented in LAMMPS. 63 In Table 1, the free energy differences calculated using our framework are compared to independent Frenkel-Ladd estimates from Prestipino et al. 64,65 The results of Prestipino et al. also include the finite size corrections of Polson et al. 66 The finite size corrections are of the order 0.01 kB T per particle, the paths to and from the Einstein crystal require on the order 100 kB T per particle, and the overall free energy difference between the polymorphs is of the order 0.1 kB T per particle. As shown in Table 1 the results predicted from our framework are similar to those from the more difficult calculations by Prestipino et al. The free energies are decomposed into entropic and energetic parts in Table 2. The free energy diabats for FCC and BCC phases of Gaussian core solids are shown in Figure 4. The raw ρa (∆E) and ρb (∆E) data is shown in Figure 5 for simulations with various values of N from 1024 to 5488 atoms. ∆E is an extensive variable, so the centers of the ρF CC (∆E)


Page 10 of 20

Page 11 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Free energy differences ∆F = FF CC − FBCC from quadratic interpolations at several densities and temperatures. The FCC and BCC phases each have N=5488 atoms. Non-harmonic interpolation gives results that deviate from quadratic interpolation results by 16% on average (See Section 1 of Supporting Information). ρσ 3 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.30 0.30

β 333.33 303.03 285.71 270.27 263.16 166.67 125.00 500.00 250.00

β∆F/N 0.18 0.17 0.16 0.15 0.15 0.09 0.06 0.32 0.15

β∆F/N 64 0.19 0.18 0.17 0.16 0.16 0.10 0.07 0.39 0.18

¯ and entropic (∆S) parts of the free energy difference Table 2: The energetic (∆E) (∆F ) for all the densities and temperatures reported in Table 1. ρσ 3 0.24 0.24 0.24 0.24 0.24 0.24 0.24 0.30 0.30

β 333.33 303.03 285.71 270.27 263.16 166.67 125.00 500.00 250.00

¯ β∆E/N 0.17 0.15 0.14 0.13 0.13 0.09 0.07 0.29 0.16

∆S/kB N -0.01 -0.02 -0.02 -0.02 -0.02 0.00 0.01 -0.03 0.01

and ρBCC (∆E) will be separated by an amount proportional to system size N . Fluctuations of ∆E within the FCC and BCC polymorphs also increase with N , but the ratio of the standard deviation to the mean decreases as N −1/2 leading to narrower diabats on a per particle basis. Therefore a small system size would seem to facilitate the interpolation. To the contrary, finite size effects are evident in the shifting peaks of ρF CC (∆E) and ρBCC (∆E) for simulations with N < 4394 atoms. The legend in Figure 5 contains the free energy



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4: Free energy diabats constructed from raw data (black, lower) for the FCC and BCC phases from an NVT simulation at ρσ 3 = 0.24 and β = 285.71. difference predicted for each size. The free energy difference approximately converges for N ≥ 4394.

Figure 5: The distributions ρF CC (∆E) and ρBCC (∆E) for values of N between 1024 atoms and 5488 atoms. The unsampled ∆E values from -0.15 to 0.30 kB T per particle are the interval in which the diabats must be interpolated.

Polymorph prediction studies rank the stabilities of polymorphs on the basis of lattice energies 24,25,28 at 0 K, i.e. with u ~ = ~0 so that entropic contributions are completely neglected. At non-zero temperatures, entropic contributions could alter the stability hierarchy from 0K calculations. This work introduced and illustrated a simple diabat interpolation method to 12 ACS Paragon Plus Environment

Page 12 of 20

Page 13 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


estimate polymorph free energy differences. The method includes anharmonic effects and entropic contributions. It requires only two unbiased simulations, so it may become possible to estimate polymorph free energy differences directly from ab initio MD simulations 67,68 when force fields are not available. The method presented in this work is limited to free energy differences at constant density, but future work will present extensions to polymorphs with different densities and direct interpolation routes to Gibbs free energy differences. We anticipate that the new method will perform well for metals, intermetallics, metal oxides, and binary salts. It is less clear whether the method will work for molecular crystals that have both stiff bonds and softer non-bonded interactions.

Supporting Information Non-parabolic interpolation results.

Acknowledements This work was supported by the National Science Foundation award 1465289 in the Division of Theoretical Chemistry. The authors thank David Quigley for stimulating discussions.

References (1) Mitscherlich, E. Sur La Relation qui existe entre la forme cristalline et les proportions chimiques. Ann. Chim. Phys. 1821, 19, 350. (2) Reddy, C. M.; Basavoju, S.; Desiraju, G. R. Sorting of polymorphs based on mechanical properties. Trimorphs of 6-chloro-2,4-dinitroaniline. Chem. Commun. (Camb.) 2005, 2439–2441. (3) Cocca, M.; Lorenzo, M. L. D.; Malinconico, M.; Frezza, V. Influence of crystal poly-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

morphism on mechanical and barrier properties of poly(L-lactic acid). Eur. Polym. J 2011, 47, 1073–1080. (4) Brinkmann, M.; Gadret, G.; Muccini, M.; Taliani, C.; Masciocchi, N.; Sironi, A. Correlation between Molecular Packing and Optical Properties in Different Crystalline Polymorphs and Amorphous Thin Films of mer -Tris(8-hydroxyquinoline)aluminum(III). J. Am. Chem. Soc 2000, 122, 5147–5157. (5) Mo, S. D.; Ching, W. Y. Electronic and optical properties of three phases of titanium dioxide: Rutile, anatase, and brookite. Phys. Rev. B 1995, 51, 13023–13032. (6) Beckmann, W.; Boisteile, R.; Sate, K. Solubility of the A , B , and C Polymorphs of Stearic Acid in Decane , Methanol , and Butanone. J. Chem. Eng. Data 1984, 29, 211–214. (7) Pudipeddi, M.; Serajuddin, A. T. M. Trends in solubility of polymorphs. J. Pharm. Sci. 2005, 94, 929–939. (8) Lehmann, O. Ueber physikalische isomerie. Zeitschrift f¨ ur Kristallographie-Crystalline Materials 1877, 1, 97–131. (9) Wille, R. L.; Lutton, E. S. Polymorphism of cocoa butter. J. Am. Oil Chem. Soc 1966, 43, 491–496. (10) Yanagisawa, K.; Ovenstone, J. Crystallization of Anatase from Amorphous Titania Using the Hydrothermal Technique: Effects of Starting Material and Temperature. J. Phys. Chem. B 1999, 103, 7781–7787. (11) Robinson, D. M.; Go, Y. B.; Mui, M.; Gardner, G.; Zhang, Z.; Mastrogiovanni, D.; Garfunkel, E.; Li, J.; Greenblatt, M.; Dismukes, G. C. Photochemical water oxidation by crystalline polymorphs of manganese oxides: Structural requirements for catalysis. J. Am. Chem. Soc 2013, 135, 3494–3501. 14 ACS Paragon Plus Environment

Page 14 of 20

Page 15 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(12) Srinivasan, K. Crystal growth of and glycine polymorphs and their polymorphic phase transformations. J. Cryst. Growth 2008, 311, 156 – 162. (13) Yu, L. Polymorphism in molecular solids: An extraordinary system of red, orange, and yellow crystals. Acc. Chem. Res. 2010, 43, 1257–1266. (14) Cruz-Cabeza, A. J.; Reutzel-Edens, S. M.; Bernstein, J. Facts and fictions about polymorphism. Chem. Soc. Rev. 2015, 44, 8619–8635. (15) Haleblian, J.; McCrone, W. Pharmaceutical applications of polymorphism. J. Pharm. Sci. 1969, 58, 911–929. (16) Chieng, N.; Rades, T.; Aaltonen, J. An overview of recent studies on the analysis of pharmaceutical polymorphs. J. Pharm. Biomed. Anal. 2011, 55, 618–644. ¨ Peterson, M. L.; Remenar, J. F.; Read, M. J.; (17) Morissette, S. L.; Almarsson, O.; Lemmo, A. V.; Ellis, S.; Cima, M. J.; Gardner, C. R. High-throughput crystallization: Polymorphs, salts, co-crystals and solvates of pharmaceutical solids. Adv. Drug Deliv. Rev. 2004, 56, 275–300. (18) Navrotsky, A. Energetic clues to pathways to biomineralization: Precursors, clusters, and nanoparticles. Proc. Natl. Acad. Sci. USA 2004, 101, 12096–12101. (19) Meldrum, F. C.; Cölfen, H. Controlling mineral morphologies and structures in biological and synthetic systems. Chem. Rev. 2008, 108, 4332–4432. (20) Grimvall, G.; Ebbsj, I. Polymorphism in Metals I. Vibrational Free Energy. Phys. Scripta 1975, 12, 168. (21) Jayaraman, A.; Klement, W.; Kennedy, G. C. Melting and Polymorphism at High Pressures in Some Group IV Elements and III-V Compounds with the Diamond/Zincblende Structure. Phys. Rev. 1963, 130, 540–547.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(22) Stack, A. G.; Gale, J. D.; Paolo, R. Restricted access Virtual Probes of Mineral-Water Interfaces: The More Flops, the Better! Elements 2013, 9, 211–216. (23) Green, H. W. Shearing instabilities accompanying high-pressure phase transformations and the mechanics of deep earthquakes. Proc. Natl. Acad. Sci. USA 2007, 104, 9133– 9138. (24) Beyer, T.; Day, G. M.; Price, S. L. The Prediction , Morphology , and Mechanical Properties of the Polymorphs of Paracetamol. J. Am. Chem. Soc 2001, 5086–5094. (25) Price, S. Computed Crystal Energy Landscapes for Understanding and Predicting Organic Crystal Structures and Polymorphism. Acc. Chem. Res. 2008, 42, 117–126. (26) Pickard, C. J.; Needs, R. J. Ab initio random structure searching. J. Phys. Condens. Matter 2011, 23, 053201. (27) Neumann, M. A.; van de Streek, J.; Fabbiani, F. P. A.; Hidber, P.; Grassmann, O. Combined crystal structure prediction and high-pressure crystallization in rational pharmaceutical polymorph screening. Nat. Commun. 2015, 6, 7793. (28) Vasileiadis, M.; Pantelides, C. C.; Adjiman, C. S. Prediction of the crystal structures of axitinib, a polymorphic pharmaceutical molecule. Chem. Eng. Sci. 2015, 121, 60–76. (29) Frenkel, D.; Ladd, A. J. C. New Monte Carlo method to compute the free energy of arbitrary solids. Application to the fcc and hcp phases of hard spheres. J. Chem. Phys. 1984, 81, 3188–3193. (30) Martonák, R.; Laio, A.; Parrinello, M. Predicting crystal structures: the ParrinelloRahman method revisited. Phys. Rev. Lett. 2003, 90, 075503. (31) Dybeck, E. C.; Schieber, N. P.; Shirts, M. R. Effects of a More Accurate Polarizable Hamiltonian on Polymorph Free Energies Computed Efficiently by Reweighting PointCharge Potentials. J. Chem. Theory Comput. 2016, 12, 3491–3505. 16 ACS Paragon Plus Environment

Page 16 of 20

Page 17 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(32) Tan, T. B.; Schultz, A. J.; Kofke, D. A. Efficient calculation of temperature dependence of solid-phase free energies by overlap sampling coupled with harmonically targeted perturbation. J. Chem. Phys. 2010, 133 . (33) Bruce, A.; Wilding, N.; Ackland, G. Free Energy of Crystalline Solids: A Lattice-Switch Monte Carlo Method. Phys. Rev. Lett. 1997, 79, 3002–3005. (34) Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.; Michaelides, A. Crystal Nucleation in Liquids: Open Questions and Future Challenges in Molecular Dynamics Simulations. Chem. Rev. 2016, 116, 7078–7116. (35) Agarwal, V.; Peters, B. Solute Precipitate Nucleation: A Review of Theory and Simulation Advances. Adv. Chem. Phys 2014, 155, 97–160. (36) Nyman, J.; Day, G. M. Static and lattice vibrational energy differences between polymorphs. CrystEngComm 2015, 17, 5154–5165. (37) Reilly, A. M.; Tkatchenko, A. Role of dispersion interactions in the polymorphism and entropic stabilization of the aspirin crystal. Phys. Rev. Lett. 2014, 113, 1–5. (38) Yu, T. Q.; Chen, P. Y.; Chen, M.; Samanta, A.; Vanden-Eijnden, E.; Tuckerman, M. Order-parameter-aided temperature-accelerated sampling for the exploration of crystal polymorphism and solid-liquid phase transitions. J. Chem. Phys. 2014, 140, 1–13. (39) Moustafa, S. G.; Schultz, A. J.; Kofke, D. A. A comparative study of methods to compute the free energy of an ordered assembly by molecular simulation. J. Chem. Phys. 2013, 139, 0–10. (40) Rickman, J. M.; LeSar, R. Free-energy Calculations in Materials Research. Annu. Rev. Mater. Res. 2002, 32, 195–217. (41) Gavezzotti, A. Ten years of experience in polymorph prediction: what next? CrystEngComm 2002, 4, 343–347. 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(42) Baroni, S.; de Gironcoli, S.; Dal Corso, A.; Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 2001, 73, 515– 562. (43) Day, G. M.; Price, S. L. A Nonempirical Anisotropic Atom-Atom Model Potential for Chlorobenzene Crystals. J. Am. Chem. Soc. 2003, 125, 16434–16443. (44) Carrier, P.; Wentzcovitch, R.; Tsuchiya, J. First-principles prediction of crystal structures at high temperatures using the quasiharmonic approximation. Phys. Rev. B 2007, 76, 1–5. (45) Pamuk, B.; Soler, J. M.; Ram´ırez, R.; Herrero, C. P.; Stephens, P. W.; Allen, P. B.; Fernández-Serra, M.-V. Anomalous Nuclear Quantum Effects in Ice. Phys. Rev. Lett. 2012, 108, 193003. (46) Heit, Y. N.; Beran, G. J. O. How important is thermal expansion for predicting molecular crystal structures and thermochemistry at finite temperatures? Acta Crystallogr. Sect. B-Struct. Sci. 2016, 72, 514–529. (47) Vega, C.; Noya, E. G. Revisiting the Frenkel-Ladd method to compute the free energy of solids: The Einstein molecule approach. J. Chem. Phys. 2007, 127 . (48) Yu, T. Q.; Tuckerman, M. E. Temperature-accelerated method for exploring polymorphism in molecular crystals based on free energy. Phys. Rev. Lett. 2011, 107, 3–6. (49) Bruce, A. D.; Jackson, A. N.; Ackland, G. J.; Wilding, N. B. Lattice-switch Monte Carlo method. Phys. Rev. E 2000, 61, 906–919. (50) Jackson, A. N.; Bruce, A. D.; Ackland, G. J. Lattice-switch Monte Carlo method: Application to soft potentials. Phys. Rev. E 2002, 65, 1–12. (51) Wilding, N. B.; Magee, J. E. Phase behavior and thermodynamic anomalies of coresoftened fluids. Phys. Rev. E 2002, 66, 1–14. 18 ACS Paragon Plus Environment

Page 18 of 20

Page 19 of 20

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(52) Bridgwater, S.; Quigley, D. Lattice-switching Monte Carlo method for crystals of flexible molecules. Phys. Rev. E 2014, 90, 1–7. (53) Underwood, T. L.; Ackland, G. J. monteswitch : A package for evaluating solid-solid free energy differences via lattice-switch Monte Carlo. arXiv preprint arXiv:1609.04329 2016, (54) Blumberger, J.; Tavernelli, I.; Klein, M. L.; Sprik, M. Diabatic free energy curves and coordination fluctuations for the aqueous Ag + Ag 2+ redox couple: A biased BornOppenheimer molecular dynamics investigation. J. Chem. Phys. 2006, 124, 0–12. (55) Hwang, J.; Warshel, A. Microscopic Examination of Free-Energy Relationships for Electron Transfer in Polar Solvents. J. Am. Chem. Soc 1987, 109, 715–720. (56) Blumberger, J.; Sprik, M. Quantum versus classical electron transfer energy as reaction coordinate for the aqueous Ru2+/Ru3+ redox reaction. Theor. Chem. Acc. 2006, 115, 113–126. (57) Marcus, R. A. Chemical and Electrochemical Electron-Transfer Theory. Annu. Rev. Phys. Chem. 1964, 15, 155. (58) Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo data. J. Comput. Phys 1976, 22, 245–268. (59) Milton, J. S.; Arnold, J. C. Introduction to probability and statistics: principles and applications for engineering and the computing sciences; McGraw-Hill, Inc., 2002. (60) Stillinger, F. H. Phase transitions in the Gaussian core system. J. Chem. Phys. 1976, 65, 3968. (61) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511–519.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(62) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695–1697. (63) Plimpton, S. Fast Parallel Algorithms for Short Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (64) Prestipino, S.; Saija, F.; Giaquinta, P. V. Phase diagram of softly repulsive systems: The Gaussian and inverse-power-law potentials. J. Chem. Phys. 2005, 123 . (65) Prestipino, S.; Saija, F.; Giaquinta, P. V. Phase diagram of the Gaussian-core model. Phys. Rev. E 2005, 71, 050102. (66) Polson, J. M.; Trizac, E.; Pronk, S.; Frenkel, D. Finite-size corrections to the free energies of crystalline solids. J. Chem. Phys. 2000, 112, 5339. (67) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: atomistic simulations of condensed matter systems. Wiley Interdiscip Rev Comput Mol Sci 2014, 4, 15–25. (68) Soler, J. M.; Artacho, E.; Gale, J. D.; Garc´ıa, A.; Junquera, J.; Ordejón, P.; SánchezPortal, D. The SIESTA method for ab initio order-N materials simulation. J. Phys. Condens. Matter 2002, 14, 2745.


Page 20 of 20

Diabat Interpolation for Polymorph Free-Energy Differences

Recommend Documents