Crystal Polymorph Search in the NPT Ensemble via a Deposition/Sublimation Alchemical Path

Abstract

The formulation of active pharmaceutical ingredients involves discovering stable crystal packing arrangements or polymorphs, each of which has distinct pharmaceutically relevant properties. Traditional experimental screening techniques utilizing various conditions are commonly supplemented with in silico crystal structure prediction (CSP) to inform the crystallization process and mitigate risk. Predictions are often based on advanced classical force fields or quantum mechanical calculations that model the crystal potential energy landscape but do not fully incorporate temperature, pressure, or solution conditions during the search procedure. This study proposes an innovative alchemical path that utilizes an advanced polarizable atomic multipole force field to predict crystal structures based on direct sampling of the NPT ensemble. The use of alchemical (i.e., nonphysical) intermediates, a novel Monte Carlo barostat, and an orthogonal space tempering bias combine to enhance the sampling efficiency of the deposition/sublimation phase transition. The proposed algorithm was applied to 2-((4-(2-(3,4-dichlorophenyl)ethyl)phenyl)amino)benzoic acid (Cambridge Crystallography Database Centre ID: XAFPAY) as a case study to showcase the algorithm. Each experimentally determined polymorph with one molecule in the asymmetric unit was successfully reproduced via approximately 1000 short 1 ns simulations per space group where each simulation was initiated from random rigid body coordinates and unit cell parameters. Utilizing two threads of a recent Intel CPU (a Xeon Gold 6330 CPU at 2.00 GHz), 1 ns of sampling using the polarizable AMOEBA force field can be acquired in 4 h (equating to more than 300 ns/day using all 112 threads/56 cores of a dual CPU node) within the Force Field X software (https://ffx.biochem.uiowa.edu). These results demonstrate a step forward in the rigorous use of the NPT ensemble during the CSP search process and open the door to future algorithms that incorporate solution conditions using continuum solvation methods.

Figure 1

Smooth alchemical path connects vacuum (λ = 0) and crystalline (λ = 1) states. During simulations in the NPT ensemble, a random walk along state variable λ facilitates sampling of both atomic coordinates and lattice parameters. Here, the path sampled for compound XXIII in the P1 space group is depicted.

Figure 2

All NPT snapshots for the two observed experimental space groups are plotted (P1 and P2₁/c). The dotted black lines denote cutoffs that prioritize high density (>1.25 g/cm³) and low potential energy (<10 kcal/mol). Panel (A) depicts the predicted structures based on minimization using AP20 parameters, whereas panel (B) was minimized with AP23 parameters. Panel (C) shows predictions after minimization with “coarse” DFT. A retrospective energy window is labeled on the right-hand side in each plot that corresponds to a relative energy that preserves the predicted polymorphs that match the experiment.

Figure 3

Shown is the density vs “precise” relative DFT-D potential energy for polymorphs discovered via the alchemical NPT pipeline and those determined experimentally. Predicted structures that match the experiment are labeled with an arrow along with the rmsd₂₀ for their respective experimental polymorph.

Figure 4

Dichlorobenzene group of compound XXIII was superposed for each saved snapshot to assess the dihedral angles sampled during the CSP procedure.

Figure 5

PAC rmsd₂₀ values for comparing simulation snapshots to experiments during the AP20 simulations that generated crystal packings close to experiment: (A) XAFPAY with space group P2₁/c, (B) XAFPAY01 with space group P1, and (C) XAFPAY03 with space group P2₁/c. The rmsd₂₀ for unminimized structures (i.e., molecular dynamics snapshots) are represented by a solid black line, whereas the rmsd₂₀ after crystal minimization is displayed by a blue dashed line. Snapshots that were eventually promoted to DFT minimization are highlighted with a blue circle.

Figure 6

Plots illustrating the OST sampling approach using the compound XXIII simulation that produced the lowest potential energy structure (based on “precise” DFT-D). Panel (A) denotes the ensemble average partial derivative of the potential energy with respect to λ (given by ⟨∂U/∂λ⟩) and its integration over the phase transition path to yield the deposition free energy difference. Panels (B,C) are contour plots of the total OST bias and only the 2D component of the bias, respectively.

Figure 7

Predicted structures that obtained the lowest relative energy based on “precise” DFT-D are compared to the AP20 relative energies. The diamonds represent the experimental polymorphs, while the circles are structures from the alchemical NPT search. Predicted structures with the lowest rmsd₂₀ to each experimental polymorph are emphasized with an ‘X’ that is the same color as the experimental polymorph.

Figure 8

Polymorphs with a relative lattice energy based on “precise” DFT-D within 2 kcal/mol of the lowest energy structure were reminimized with the AMOEBA Poltype2 2023 parameters under three polarization models: (A) no polarization. (B) Direct polarization. (C) Mutual polarization.