A Cannibalistic Approach to Grand Canonical Crystal Growth - Journal

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

A Cannibalistic Approach to Grand Canonical Crystal Growth Tarak Karmakar,†,§ Pablo M. Piaggi,‡,§ Claudio Perego,#,†,§ and Michele Parrinello*,†,§ †

Department of Chemistry and Applied Biosciences, ETH Zurich, c/o USI Campus, Via Giuseppe Buffi 13, CH-6900, Lugano, Ticino, Switzerland ‡ Theory and Simulation of Materials (THEOS), École Polytechnique Fiédérale de Lausanne, CH-1015 Lausanne, Switzerland § Facoltià di Informatica, Instituto di Scienze Computationali, Università della Svizzera Italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Ticino, Switzerland # Polymer Theory Department, Max Planck Institute for Polymer Research, Ackermannweg 10, 55128 Mainz, Germany S Supporting Information *

ABSTRACT: Canonical molecular dynamics simulations of crystal growth from solution suffer from severe finite-size effects. As the crystal grows, the solute molecules are drawn from the solution to the crystal, leading to a continuous drop in the solution concentration. This is in contrast to experiments in which the crystal grows at an approximately constant supersaturation of a bulk solution. Recently, Perego et al. [J. Chem. Phys. 2015, 142, 144113] showed that in a periodic setup in which the crystal is represented as a slab, the concentration in the vicinity of the two surfaces can be kept constant while the molecules are drawn from a part of the solution that acts as a molecular reservoir. This method is quite effective in studying crystallization under controlled supersaturation conditions. However, once the reservoir is depleted, the constant supersaturation conditions cannot be maintained. We propose a variant of this method to tackle this depletion problem by simultaneously dissolving one side of the crystal while letting the other side grow. A continuous supply of particles to the solution due to the crystal dissolution maintains a steady solution concentration and avoids reservoir depletion. In this way, a constant supersaturation condition can be maintained for as long as necessary. We have applied this method to study the growth and dissolution of urea crystal from water solution under constant supersaturation and undersaturation conditions, respectively. The computed growth and dissolution rates are in good agreement with those obtained in previous studies.

1. INTRODUCTION Understanding crystallization is of paramount interest in the field of material, pharmaceutical, and environmental sciences.1,2 Over the years, computer simulations have emerged as effective tools in the study of crystal growth and predicting polymorphism.3,4 Thus, molecular dynamics (MD) and Monte Carlo (MC) techniques are routinely employed to obtain crystallization atomic details that are often difficult to get from experiments.5−7 In an experiment, crystallization is usually carried out at a constant solution concentration. To mimic the crystallization process on the computer, one usually considers model systems containing a few thousand molecules. Such small systems suffer from severe finite-size effects.8,9 One limitation that concerns us here is the fact that since most MD simulations are carried out in the canonical ensemble, as the crystal grows, the solution concentration decreases. This dramatically affects the crystal growth and leads to an erroneous estimation of the growth rate, unless proper corrections are introduced.9−13 Numerical methods such as Grand-Canonical (GC) MC and MD have been developed to tackle this issue.14−17 These methods are based on particle insertion and deletion steps. However, in a dense fluid, these methods become inefficient © XXXX American Chemical Society

since the probability of inserting a molecule is very low. The recently developed adaptive-resolution MD technique provides a solution to this problem by allowing the exchange of particles between a low-resolution coarse-grained and a high-resolution atomistic region.18 However, coupling of the different resolution regions is a challenging task,19 which can impose limitations. Despite several attempts to solve the finite-size problem, the search for a simpler approach that could circumvent this issue is still going on. Very recently, Perego et al.20 designed a constant chemical potential molecular dynamics (CμMD) method in which the concentration of the solution in the immediate surroundings of the growing crystal is kept constant by means of an applied external force, while the rest of the solution acts as a molecular reservoir. This leads to a steady crystal growth under a constant chemical potential. This method has been applied to study the growth of a urea crystal from its supersaturated aqueous solution. The growth rates of the (001) and (110) surfaces of the urea crystal have been calculated, and their values agree quite well with those reported Received: February 21, 2018 Published: March 28, 2018 A

DOI: 10.1021/acs.jctc.8b00191 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 1. A schematic description of the simulation box. The crystal is placed in the middle. Partitions are shown in different colors and marked on top. Urea molecules are depicted in blue sticks, while water molecules are shown in red spheres. Abbreviations and symbols introduced in the figure are discussed in the text.

in previous studies.21 Although this method has opened up a new way of studying crystallization under controlled supersaturation conditions,22 it suffers from a significant limitation. As the crystal grows, more urea molecules are drawn from the reservoir, leading to its depletion. When the reservoir is strongly depleted, the effectiveness of the CμMD scheme comes to a halt. We have developed a variant of the CμMD method to overcome the depletion problem and maintain a constant solution supersaturation near the crystal for a theoretically indefinite time. In this method, the crystal grows on one side in the presence of a supersaturated solution and dissolves on the other side that is exposed to an undersaturated solution. During the crystal growth, the solute molecules are drawn from the solution and deposited on the crystal surface, while the dissolution on the other side of the crystal slab leads to a constant supply of solute particles to the reservoir, thus avoiding its depletion. Our work has been inspired by ref 23, in which a constant concentration gradient has been set up across a membrane. We have applied the method presented in this work to study the growth of the urea (001) surface from a supersaturated water solution. We are able to maintain a steady state for all the simulation duration. This allows the growth of several crystal layers to be simulated for a long time, thus improving statistics. Furthermore, the dissolution kinetics of urea under constant undersaturation (sink condition) can be investigated at the same time. An excellent agreement between growth and dissolution rates of the urea crystal with those reported earlier in the literature21,24 demonstrates the usefulness of our method.

ct =

N

1 V CR

∑ θ (z j ) (1)

j=1

where VCR is the volume of the CR and ⎧1, if zj ∈ CR θ (z j ) = ⎨ ⎪ ⎩ 0, otherwise ⎪

(2)

is an indicator function that identifies the particles belonging to a CR. In practice, θ(z) is switched continuously from 1 to 0 at the CR boundaries to avoid sudden jumps in density.20 The external force Fc(z) has the following expression: F c(z) = k(ct − c)Gσ (z , ZF )

(3)

where k is the force constant, c is the target concentration, while ct indicates its instantaneous value, and Gσ(z, ZF) is a bellshaped function that localizes the action of Fc(z) around the position ZF = zI + DF. The force, Fc, acts as a harmonic-like restraint which maintains the solution concentration (ct) close to a targeted value (c) in the CRs, while the region beyond DF is used as a molecular reservoir. Gσ(z, ZF) in eq 3 is defined as follows: Gσ (z − ZF ) =

−1 ⎛ z − Z F ⎞⎤ 1 ⎡ ⎟⎥ ⎢1 + cosh⎜⎝ 4σ ⎣ σ ⎠⎦

(4)

Gσ is nonzero in the vicinity of the force center (z ∼ ZF) with an intensity peak proportional to σ−1 and a width proportional to σ. Depending on the force and solution features, the action of this restraint will induce a supersaturation gradient across ZF, which extends for a width σF (force region, FR shown in red color, Figure 1). Note that σ is different from σF. The CRs are defined in such a way that they do not overlap with the FR. Analogously, the presence of the crystal surface determines a supersaturation gradient within the solution that extends to a distance ξ from the interface (transition region, TR shown in green color, Figure 1). As shown in Figure 1, the CRs are defined so that they do not overlap with the TR as well. With this setting, the concentration within the CR has a nearly flat profile, thus mimicking a bulk solution. The method has been implemented in the following way. The first step is to locate the two crystal−solution interfaces (zI). In order to do that, the simulation box is divided into bins of size δz = 0.075 nm, and the solvent density in each bin is calculated. The solvent density is zero inside the crystal, and we identify the crystal−solution interfaces where the solvent density exceeds the threshold value of 5 nm−3. Once the

2. COMPUTATIONAL METHODS Constant Chemical Potential Molecular Dynamics (CμMD). In this section, we briefly discuss the principles underlying the CμMD method. As discussed in ref 20, in this method, the simulation box is partitioned into different regions (Figure 1). Periodic boundary conditions are imposed in all three directions; z is the direction perpendicular to the growing crystal surface. For clarity, we place the crystal slab in the middle of the simulation box. At each crystal side, a control region (CR, shown in blue color in the schematic, Figure 1) is defined, in which the solution concentration (c) is kept constant by an external force (Fc) applied at a distance DF from the interface, zI (as shown in Figure 1). The instantaneous solution concentration ct in the CR is evaluated as B


Article

Journal of Chemical Theory and Computation Table 1. Summary of the Simulated Systemsa

interfaces are found, the CR boundaries and the position of the external force ZF need to be defined. In order to do that, we use the values of ξ and σF that have been estimated based on a previous study20 and also from new test simulations. As the crystal grows and dissolves, the CRs and ZF are maintained at a constant distance from the interfaces, the position of which (zI) is updated on-the-fly. In the original implementation of CμMD, the forces impose the same solution supersaturation at both crystal sides, i.e., maintaining a unique c in the two CRs. We instead impose a gradient in the chemical potential between the two CRs, by keeping a supersaturated solution in one CR (CRS), with target concentration cS, and an undersaturated solution in the other CR (CRU), with a target concentration, cU.

type Unbiased (320 K)

Biased (300 K)

systems

cS (nm−3)

cU (nm−3)

A1 A2 A3 B1 B2 B3

1.5 2.0 2.4 1.5 2.0 2.4

0.4 0.2 0.2 0.5 0.5 0.5

Simulations carried out at 320 and 300 K are referred to as “A” and “B,” respectively. The cS and cU are supersaturation and undersaturation of the solution present in CRS and CRU, respectively.

a

longer obtain information about the dissolution rate. At 300 K, we performed three simulations (B1, B2, and B3) under the same supersaturation conditions as in the 320 K case (Table 1).

3. A CASE STUDY: UREA CRYSTAL GROWTH System Preparation. The method has been applied to study the growth of urea crystal in its fastest growing (001) surface. We considered a system of 679 urea and 1757 water molecules in a box of dimension 2.68 × 2.68 × 13.8 nm3. Out of the 679 urea molecules, 250 were arranged to form a crystal slab of five layers. The rest were dissolved in the solution (Figure 1). The initial configuration in all simulations was prepared as follows. First, a configuration was extracted from ref 20, then a concentration gradient was applied to reach the desired high (cS) and low (cU) solution concentration in CRS and CRU, respectively. Simulation Details. All simulations were carried out using the Gromacs-5.1.4 code25−28 patched with a private version of Plumed2.29 The simulations were run in the NVT ensemble. The stochastic velocity rescaling (Bussi−Donadio−Parrinello) thermostat30 was used to control the temperature. A time step of 2 fs was used to integrate the equations of motion. The covalent bonds involving hydrogen atoms were constrained using the LINCS algorithm.31 A cutoff of 1 nm was used for both van der Waals and short-range Coulomb interactions. Long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) method. Amber force-field parameters were used to model urea along with the TIP3P model for water molecules.32,33 The visual molecular dynamics (VMD)34 software was used to visualize the trajectories and prepare a movie. We describe here two types of simulations. At 320 K, we started imposing a supersaturation concentration (cS) of 1.5 nm−3 (system A1). After some experimentation, we found that in the presence of an undersaturated solution of (cU) of 0.4 nm−3, the system was able to remain in a steady state in which the rate of growth compensates that of dissolution. We performed two additional simulations (A2 and A3) of this type at cS = 2.0 nm−3 and 2.4 nm−3. In these cases, the cU was maintained at a much lower value of 0.2 nm−3 (Table 1) such that the enhanced growth rate was compensated by the increased dissolution rate. When the temperature was brought down to 300 K, the growth rate increased, and dissolution at the same time decreased (the results obtained from this test simulation are provided in section 1 of SI). Thus, it was not possible to find undersaturation conditions low enough for the dissolution to compensate the growth (Figure S1). In order to accelerate the dissolution, we applied an extra bias potential (Vb) to the dissolving layer (Figure S2, see SI for details). In this way, we were able to maintain a steady growth for the entire simulation time. Of course, since we had added the bias, we could no

4. RESULTS AND DISCUSSION Unbiased Simulations at 320 K. We first describe the results obtained from the unbiased simulations performed at a temperature of 320 K. As the dynamics progressed, the crystal has grown on the side exposed to the supersaturated solution and dissolved on the side that was facing the undersaturation solution. To keep track of the number of crystalline urea molecules, we used the crystal fingerprint variable (Ncr)4,35 that distinguishes the crystalline molecules and those present in the solution. In the A1 and A2 simulation trajectories, the crystal keeps on average a constant number of molecules during the entire simulation (Figure S3a). When we increased the cS further to 2.4 nm−3, the steady state condition was not perfect, since the growth was faster than the dissolution. However, as discussed in the following, the mismatch between the growth and dissolution was small enough to maintain the constant supersaturation conditions over a long enough time to measure the growth rate. We shift our attention to the solution concentration in the CRs. Both ctS and ctU fluctuate near their corresponding target concentrations cS and cU, respectively, demonstrating the effectiveness of the method in maintaining a constant CR concentration for the whole simulation duration (Figure S4). Unlike the A1 and A2 trajectories, in A3, a faster growth rate compared to that of the dissolution results in a small deviation from the perfect steady state. These deviations however are so small that their effect cannot be seen in the concentrations of the CRs (Figure S4). We then calculated the growth rates from the simulation trajectories. At first, an imaginary plane was placed at a distance of 1.0 nm from the crystal−solution interface (zI), at the supersaturated side of the box. The number of urea molecules that cross the plane from left to right (N+) and vice versa (N−) were counted. The total number of urea molecules (ΔN = N+ − N−) that eventually crossed the plane and subsequently got deposited on the growing crystal surface at a given time was calculated from ΔN vg = nstepsΔt (5) where vg is the growth speed, nsteps is the number of steps, and Δt is the integration time step. The growth rate was calculated from the slope of the ΔN vs time plot (Figure 2). The statistical errors in the rates were estimated using the block-average technique.36 The ability to reach this steady state is determined by the balance between the increase of the dissolution rate and C


Article


dissolution rates obtained from the linear fit of the ΔN curves at cU = 0.4 nm−3 and 0.2 nm−3 are found to be 1.26 ± 0.32 ns−1 and 1.50 ± 0.37 ns−1, respectively. It is interesting to note that the crystal dissolution takes place in a stepwise fashion (Figure 3), the first step of which is the detachment of a few molecules from the crystal surface resulting in the formation of defect sites followed by the entire layer dissolution, similar to the observations by Piana and Gale.37 In such a case, a linear fit could be a crude approximation. Thus, we performed an additional analysis collecting the time intervals of the dissolution of consecutive layers from four independent simulation trajectories (discussed in SI, Figure S6a). A cumulative time distribution (Figure S6b and c) fitted to an exponential model provides another estimate of the dissolution time of ∼30 ns. This corresponds to a dissolution rate of 1.67 ns−1, compatible with the values obtained from the linear fitting (Table 2). Similar analyses were not performed for the higher cU due to insufficient statistics. Biased Simulations at 300 K. Next, we present the results obtained from the biased simulations. At this lower temperature, spontaneous dissolution cannot compensate the increase in growth rate. For this reason, we help dissolution by imposing a harmonic potential (Vb) that speeds up the process. We found that it was possible to reach a steady state using the following set of spring contants (eq 5 of SI) kd = 1.0, 2.0, and 2.75 kJ/mol for B1, B2, and B3, respectively (Figure S7). Lower values of kd reflect the lower dissolution rates. Within the 100 ns simulation time (Figure S3b), several crystal layers have grown under steady conditions (Figure S8), allowing us to gather enough statistics. The growth rates obtained from the biased simulations are reported in Table 2. Within error bars, the growth rates are equal to those reported by Perego et al.20 (Figure 4).

Figure 2. Growth profiles obtained from the (a) unbiased (320 K) and (b) biased (300 K) simulation trajectories. ΔN is defined in the text.

decrease of the growth rate. With decreasing the temperature, the growth rate increases, while the dissolution becomes slower (Table 2). Table 2. Growth and Dissolution Rates Obtained from the Unbiased (320 K) Simulations (A1, A2, and A3)a type

system

unbiased (320 K)

A1 A2 A3 B1 B2 B3

biased (300 K)

growth rate (ns−1)b

dissolution rate (ns−1)c

± ± ± ± ± ±

1.26 ± 0.32 1.50 ± 0.37 1.50 ± 0.41

1.46 2.45 3.63 2.12 3.26 5.65

0.45 0.53 0.61 0.73 0.70 0.68

a

In the case of the biased (300 K) simulations (B1, B2, and B3), only growth rates were calculated. bIt is to be noted that in the simulations carried out by Perego et al., the urea crystal grows on both sides, while in our case, it grows on one side. Thus, the reported growth rate in ref 20 is almost twice the value we have obtained in this work. cThe dissolution rates were not calculated from the biased simulation trajectories.

We followed a similar protocol to calculate the dissolution rates. The dissolution profiles are shown in Figure 3. The Figure 4. Growth rates plotted as a function of the supersaturation in CRS obtained from the unbiased (blue) and biased (red) simulation trajectories. The filled circles (along with error bars) indicate the absolute growth rates, while the lines depict the linear fits. For the sake of comparison, the growth rates reported in ref 20 are shown by the asterisk symbols.

In Figure 4, the growth rates obtained from both the unbiased and biased simulation trajectories are plotted as a function of supersaturation. A linear relationship between the growth rate and the solution supersaturation, which is generally associated with the expected rough growth mechanism,38 has been observed for the urea (001) surface growth. Of course, since we have biased the dissolving layer, no rates could be evaluated for this process. However, also in this case, the

Figure 3. Dissolution profiles obtained from the unbiased (320 K) simulation trajectories. The negative y axis indicates the loss of the crystalline molecules. Both growth and dissolution profiles are provided in Figure S5. D



■

ACKNOWLEDGMENTS The authors would like to thank Dr. Omar Valsson and Roberto Menichetti for carefully reading the manuscript and providing useful suggestions. The computational resources are provided by the CSCS, Swiss National Supercomputing Centre. The research was supported by the European Union Grant No. ERC-2014-AdG-670227/VARMET. We also acknowledge the NCCR MARVEL, funded by the Swiss National Science Foundation.

dissolution process appeared to be triggered by the creation of defects.

5. CONCLUSIONS In summary, we have developed a method to simulate crystallization from solution that is capable of maintaining a constant supersaturation in the vicinity of the crystal surface, leading to a steady growth process unaffected from solution depletion issues. The method is a variant of the CμMD scheme, in which the finite reservoir effect is removed by enforcing dissolution on one crystal side while letting the other side grow. As a result, unlike previous CμMD simulations, in the present implementation, the constant solution supersaturation is maintained for a much longer time, and in theory, it could be maintained indefinitely. The growth rates of the urea (001) surface calculated from the biased simulation trajectories at different supersaturations are in excellent agreement with that obtained from the previous CμMD simulations. Furthermore, both growth and dissolution rates are obtained from the unbiased simulation trajectories at a temperature of 320 K. The dissolution rates calculated from the unbiased simulation trajectories are comparable to those recently reported.39 In this calculation, the ability, under appropriate circumstances, to calculate dissolution rates is a welcome byproduct. Wherever a spontaneous dissolution is not able to balance the growth rate, this is not possible since dissolution needs to be accelerated by external bias. If, however, as in many pharmaceutical applications, one were interested mostly in dissolution, one could set up a simulation in which a bias potential is applied to the growing surface. Metadynamics40,41 or the recently developed variationally enhanced sampling42,43 techniques can be coupled with the proposed CμMD scheme to investigate a slow growing (and or a dissolving) crystal with significantly reduced computational cost. In our future work, we plan to test this method on such systems and further extend our investigations on more complex multicomponent systems.

■

■

REFERENCES

(1) Greiner, M.; Elts, E.; Briesen, H. Insights into pharmaceutical nanocrystal dissolution: a molecular dynamics simulation study on aspirin. Mol. Pharmaceutics 2014, 11, 3009−3016. (2) Tan, W.; Yang, X.; Duan, X.; Zhang, X.; Qian, G.; Zhou, X. Understanding supersaturation-dependent crystal growth of L-alanine in aqueous solution. Cryst. Res. Technol. 2016, 51, 23−29. (3) Anwar, J.; Boateng, P. K. Computer simulation of crystallization from solution. J. Am. Chem. Soc. 1998, 120, 9600−9604. (4) Salvalaglio, M.; Perego, C.; Giberti, F.; Mazzotti, M.; Parrinello, M. Molecular-dynamics simulations of urea nucleation from aqueous solution. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, E6−E14. (5) Piana, S.; Reyhani, M.; Gale, J. D. Simulating micrometre-scale crystal growth from solution. Nature 2005, 438, 70. (6) Reilly, A. M.; Briesen, H. Modeling crystal growth from solution with molecular dynamics simulations: approaches to transition rate constants. J. Chem. Phys. 2012, 136, 034704. (7) Reilly, A. M.; Briesen, H. A detailed kinetic Monte Carlo study of growth from solution using MD-derived rate constants. J. Cryst. Growth 2012, 354, 34−43. (8) Wedekind, J.; Reguera, D.; Strey, R. Finite-size effects in simulations of nucleation. J. Chem. Phys. 2006, 125, 214505. (9) Salvalaglio, M.; Tiwary, P.; Maggioni, G. M.; Mazzotti, M.; Parrinello, M. Overcoming time scale and finite size limitations to compute nucleation rates from small scale well tempered metadynamics simulations. J. Chem. Phys. 2016, 145, 211925. (10) Reguera, D.; Bowles, R.; Djikaev, Y.; Reiss, H. Phase transitions in systems small enough to be clusters. J. Chem. Phys. 2003, 118, 340− 353. (11) Grossier, R.; Veesler, S. Reaching one single and stable critical cluster through finite-sized systems. Cryst. Growth Des. 2009, 9, 1917− 1922. (12) Schmelzer, J. W.; Abyzov, A. S. Comments on the thermodynamic analysis of nucleation in confined space. J. NonCryst. Solids 2014, 384, 2−7. (13) Agarwal, V.; Peters, B. Nucleation near the eutectic point in a Potts-lattice gas model. J. Chem. Phys. 2014, 140, 084111. (14) Yau, D. H.; Liem, S. Y.; Chan, K.-Y. A contact cavity-biased method for grand canonical Monte Carlo simulations. J. Chem. Phys. 1994, 101, 7918−7924. (15) Papadopoulou, A.; Becker, E. D.; Lupkowski, M.; van Swol, F. Molecular dynamics and Monte Carlo simulations in the grand canonical ensemble: Local versus global control. J. Chem. Phys. 1993, 98, 4897−4908. (16) Lo, C.; Palmer, B. Alternative Hamiltonian for molecular dynamics simulations in the grand canonical ensemble. J. Chem. Phys. 1995, 102, 925−931. (17) Lynch, G. C.; Pettitt, B. M. Grand canonical ensemble molecular dynamics simulations: Reformulation of extended system dynamics approaches. J. Chem. Phys. 1997, 107, 8594−8610. (18) Potestio, R.; Fritsch, S.; Espanol, P.; Delgado-Buscalioni, R.; Kremer, K.; Everaers, R.; Donadio, D. Hamiltonian adaptive resolution simulation for molecular liquids. Phys. Rev. Lett. 2013, 110, 108301. (19) Kreis, K.; Fogarty, A. C.; Kremer, K.; Potestio, R. Advantages and challenges in coupling an ideal gas to atomistic models in adaptive resolution simulations. Eur. Phys. J.: Spec. Top. 2015, 224, 2289−2304.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00191. Results of the 300 K unbiased simulation, details of the biased potential, ΔNcr vs time plots, ctS and ctU, growth and dissolution profiles, cumulative distribution plots obtained from the A-type simulations, growth and dissolution profiles, and ctS and ctU plots of the B-type simulation trajectories (PDF) A movie demonstrating the growth and dissolution processes (MPG)

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Tarak Karmakar: 0000-0002-8721-6247 Claudio Perego: 0000-0001-8885-3080 Notes

The authors declare no competing financial interest. E


Article


(42) Valsson, O.; Parrinello, M. Variational approach to enhanced sampling and free energy calculations. Phys. Rev. Lett. 2014, 113, 090601. (43) Valsson, O.; Parrinello, M. Well-tempered variational approach to enhanced sampling. J. Chem. Theory Comput. 2015, 11, 1996−2002.

(20) Perego, C.; Salvalaglio, M.; Parrinello, M. Molecular dynamics simulations of solutions at constant chemical potential. J. Chem. Phys. 2015, 142, 144113. (21) Salvalaglio, M.; Vetter, T.; Mazzotti, M.; Parrinello, M. Controlling and predicting crystal shapes: The case of urea. Angew. Chem., Int. Ed. 2013, 52, 13369−13372. (22) Radu, M.; Kremer, K. Enhanced Crystal Growth in Binary Lennard-Jones Mixtures. Phys. Rev. Lett. 2017, 118, 055702. (23) Ozcan, A.; Perego, C.; Salvalaglio, M.; Parrinello, M.; Yazaydin, O. Concentration gradient driven molecular dynamics: a new method for simulations of membrane permeation and separation. Chem. Sci. 2017, 8, 3858−3865. (24) Salvalaglio, M.; Vetter, T.; Giberti, F.; Mazzotti, M.; Parrinello, M. Uncovering molecular details of urea crystal growth in the presence of additives. J. Am. Chem. Soc. 2012, 134, 17221−17233. (25) Berendsen, H. J.; van der Spoel, D.; van Drunen, R. GROMACS: a message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 1995, 91, 43−56. (26) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. GROMACS: fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (27) Hess, B.; Kutzner, C.; Van Der Spoel, D.; Lindahl, E. GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (28) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845. (29) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604−613. (30) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (31) Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463−1472. (32) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (33) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (34) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (35) Giberti, F.; Salvalaglio, M.; Mazzotti, M.; Parrinello, M. Insight into the nucleation of urea crystals from the melt. Chem. Eng. Sci. 2015, 121, 51−59. (36) Petersen, H.; Flyvbjerg, H. Error estimates in molecular dynamics simulations. J. Chem. Phys. 1989, 91, 461−466. (37) Piana, S.; Gale, J. D. Understanding the barriers to crystal growth: dynamical simulation of the dissolution and growth of urea from aqueous solution. J. Am. Chem. Soc. 2005, 127, 1975−1982. (38) Lovette, M. A.; Browning, A. R.; Griffin, D. W.; Sizemore, J. P.; Snyder, R. C.; Doherty, M. F. Crystal shape engineering. Ind. Eng. Chem. Res. 2008, 47, 9812−9833. (39) Anand, A.; Patey, G. N. The Mechanism of Urea Crystal Dissolution in Water from Molecular Dynamics Simulation. J. Phys. Chem. B 2018, 122, 1213. (40) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (41) Barducci, A.; Bussi, G.; Parrinello, M. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett. 2008, 100, 020603. F


A Cannibalistic Approach to Grand Canonical Crystal Growth - Journal

Recommend Documents