in Material erials s Science ce Mark E. Tuckerman Dept. of - - PowerPoint PPT Presentation
Enhanced En nced Free En Energy Based Structure ture Predicti iction on in Material erials s Science ce Mark E. Tuckerman Dept. of Chemistry and Courant Institute of Mathematical Sciences New York University, 100 Washington Square East,
Mark E. Tuckerman
New York University, 100 Washington Square East, NY 10003 NYU-ECNU Center for Computational Chemistry at NYU Shanghai 200062, China 纽约大学 - 华东师范大学计算化学中心在上海纽约大学
NEW YORK UNIVERSITY MRSEC
Crystallization and the Pharmaceutical Industry: Materials Science? Pharmaceutical Materials Science: An Active New Frontier in Materials Research
Guest Editors MRS Bulletin 31, Issue 11 (2006)
Pharmaceutical Industry
small molecule crystals “One of the continuing scandals in the physical sciences is that it remains impossible to predict the structure of even the simplest crystalline solids from their chemical composition.”
Targets of the sixth blind structure prediction test
Empirical potentials fit to SAPT(DFT) calculations
( ) ( ) (1) ( ) 1 bonded 6,8,10
ij ij
n r i j ij l l n ij ij ij ij n ij ij n ij l n ij ij
( , ) 1 !
m n r n m
r f r e m
Machine-learning (‘mathematical’) potentials Standard force fields
bonded 12 6
i j ij ij ij ij ij ij
Kohn-Sham density functional theory (plus dispersion):
2 { }
1 1 ( ) ( ) ( ) [ ] [ , ] ( ) 2 2 '
i i xc ext NN i
n n U d d E n E n U
r r R r r R R r r
2
i i j ij i
P212121
Lowest Energy Structures
P21/n P21/n Pbca C2/c Pbca
Optimize monomer geometry using DFT [PBE0 + D3, 6-311G**(d,p) Generate Unit Cells (UPACK)
25,000 generated structures 4,042 unique after clustering 107 within 10 kJ/mol of lowest energy
Run MD at exptl. P,T
Summary of 50 predicted structures based on MD averaged energies
I: P21/n II: Pbca, E = 3.46 kJ/mol III: P21/c, E = 5.37 kJ/mol IV: P21/n, E = 5.99 kJ/mol
V: P21/c, E = 6.31 kJ/mol
At 150 K: RMSD20 from experimental structure = 0.15 Å. Rank = 4 on our list Final ranking based on POLY1 SAPT(DFT) potential At 300 K: RMSD20 from experimental structure = 0.18 Å. Rank = 4 on our list Attempted predictions = 21, Times predicted = 12, Times ranked 5 or lower = 7
HIV-1 protease Polymorphism refers to the ability
crystal structures.
Ritonavir: HIV protease inhibitor
capsules, no refrigeration required
polymorph (form I to form II) on the shelf
gel caps (2002).
Ranitidine hydrochloride (Zantac) C13H22N4O3S
Bond et al. Amer. Pharma. Rev. (2007)
Acetylsalicylic acid
Rough energy landscapes in proteins and crystals
Dunitz and Scheraga PNAS 101, 14309 (2004)
Sampling via collective variables in molecular simulation
Peptides: Crystals:
Cell shape:
x x x y y y z z z
2 1
4 ( ) ( , ) (2 1)
b
N l l b lm b b m l b
Q f r Y l N
( , ) 1 1
n H n
p r
Adiabatic free energy (temperature-accelerated molecular) dynamics (AFED/TAMD)
Margliano and Vanden-Eijnden, Chem. Phys. Lett. 426, 168 (2006); J. B. Abrams and MET, J. Phys. Chem. B 112, 14752 (2008) For prediction of crystal structures: Yu and MET Phys. Rev. Lett. 107, 015701 (2011); Yu et al. J. Chem. Phys. 150, 214109 (2014)
1/2 2 { }
Write δ-functions as product of Gaussians: Canonical probability distribution and free energy surface:
2 2 1 1
n n s s
1
N
Suppose n collective variables characterize a free energy landscape of interest
1 1
n B n
Dynamically sample Gaussian centers s1,…,sn by introducing a kinetic energy for them And extending the phase space, which gives the following Hamiltonian:
Adiabatically decoupled equations of motion:
i i i i
s
s i
Under adiabatic conditions, we generate a distribution
({ }) adb 1
n s
({ }) adb 1 { } 1 1 { }
lim ln ( ,..., , , ) lim ln ( ,..., , ) ( ,..., , )
s n s s n n s
kT P s s T T T kT P s s T A s s T T
/ ({ }) adb 1 1 { }
s
T T n s n
Adiabatic free energy (temperature-accelerated molecular) dynamics (AFED/TAMD)
Other developments: 1. Combining with metadynamics used as a bias:
s
Constant pressure version with cell matrix as CVs
II(01)
X-ray powder diffraction patterns Yonetani et al. (2001).
II(98) 98 01
MURI/ARO
Free energies of benzene polymorphs Six-dimensional free energy surface analyzed using clustering based on a Gaussian mixture model.
II(98) II(01)
Cmca P41212 P21/c C2/c P21/c
(see, also, J. Yang et al. Science 345, 640 (2014))
depending on the relative orientations of the OH bonds, both having space group Pna21
1.
α → δ at pressures above 3 GPa2. Structure unknown.
from β phase above 5.6 GPa. This phase is also thought to be disordered based on Raman scattering2.
conformation recently found3. Its exact structure is unknown4.
resorcinol.5
1Robertson and Ubbelohde Proc. Royal Soc. London 167 122 (1938) 2 Rao et al. Phys. Rev. B 65 054108 (2002). 3 Kahr et al. Cryst. Growth Design 11 2070 (2011). 4 Zhu et al. Angew. Chem. Intl. Ed. (submitted) 5 J. Phys. Chem. B (2015); J. Mol. Struct. (2015).
polymorph
1
Preliminary free energy estimates for resorcinol
form δ form a/A b/A c/A a/A b/A c/A Cell length 10.53 9.53 5.60 8.20 11.65 4.88 Cell angle 90 90 90 90 90 90 Space GP Pna21 Pna21
Supercell length A
DADP TCTNB
β ( iodo- ) α ( chloro- ) ( bromo- )
Experimental Structures
Polarizable system without pre-determined force field parameters = opportunity for ab initio enhanced sampling!
1 3
CP2K 2.5: PBE+D3/DZVP-MOLOPT-(SR)-GTH
For optimized unit cells within PBE+D3, DZVP, 450 Ry in CP2K Cl α < Cl β ΔE ≈ 1.6 kJ/mol Br α > Br β ΔE ≈ 4.9 kJ/mol I α >> I β ΔE ≈ 15.8 kJ/mol
β ( iodo- ) α ( chloro- ) ( bromo- )
kinetic product – aggregation promotes co-crystal growth?
difference between the two forms for each co-crystal?
Start with nuclei Compute ,
i n
i n
Propagate nuclei a short time Δt with F Add electrons Add electrons
2
( ) 2 (0) ( ) (0)
I I I I I
t t t M R R R F
e.g. Verlet:
{ }
I I I
11.0 picosecond of ab initio crystal-AFED with CP2K for TCTNB/DADP co-crystal
Experimental structure: α New Polymorph Alternate polymorph: β
3
3
3
MURI/ARO
Crystal structures of Xenon at 2700 K and 25 GPa Using Q4 and Q6 only as collective variables Ts = 1.5 x 105 K ms = 108 amu Using Q4, Q6, and the cell matrix hd as collective variables Ts = Th = 105 K, ms = 108 amu
Experimental evidence mixed fcc-hcp structures in pressure-induced Martensitic transitions from fcc to hcp as intermediates: Jephcoat et al. Phys. Rev. Lett. 59, 2670 (1987) Cynn et al. Phys. Rev. Lett. 86, 4552 (2001)
Classical nucleation theory:
3 2
s
3 2
s s
r*
Collective variables: hd, Q4, Q6, T = Tm, Ts=Th = 107 K Potential model: Embedded Atom Model (EAM) of Mishin et al.
Free energy surfaces
One-dimensional profiles constructed using string method [Maragliano et al. JCP (2006)]
Navigating on a high-dimensional surface From a minimum, we ask for the nearest accessible saddle points. These can be found using gentlest ascent dynamics (GAD) on the FES.
A. Samanta, M. Chen, T. –Q. Yu, W. E, MET J. Chem. Phys. 140, 164109 (2014).
Activation: Combine with relaxation to next minimum:
Termed “STochastic Activation-Relaxation Technique” or START
Example: Alanine tripeptide Landmarks determined using density-based clustering methods. Minima and saddles are assigned as vertices of different colors. An edge connects a saddle vertex to a minimum vertex if a relaxation initiated at the saddle ends up at the corresponding minimum. Graph representation of free energy surfaces
Full network diagram for met-enkephalin using 10 Ramachandran angles as CVs: 1081 minima, 1431 saddles Graph nodes sorted based on Sketch Map analysis
2 2
ij ij i j
If we want to compute free energies, use the minima obtained in a START run to indicate where a set of bins should be placed and use AFED sampling to obtain the populations in those bins. Can be adapted for crystal structure prediction! [M. Chen, T. –Q. Yu, MET (in prep)] Can handle explicit solvent!
Experimental structure first found in 1 ns using orientational CVs. RMSD4 of carbons = 0.16 Å, full RMSD4 = 0.44 Å using a simple, universal force field
MET (in prep)
1. Protocol for crystal structure prediction including specialized force field generation [SAPT(DFT)] combined with molecular dynamics refinement. 2. Enhanced sampling approaches further allow for polymorph prediction and ranking based on free energy, which is the proper figure of merit. 3. Surface navigation combined with enhanced sampling looks to be a promising approach going forward, however, it is important to have a robust set of internal collective variables. 4. Considerable challenges remain including:
5. From Nature press release based on CCDC 6th blind test, “It is fair to claim that, to a large extent, this blind test has shown that the problem of organic structure prediction has been solved.” – M. Neumann. This claim might be exaggerated, and the field of crystal structure prediction far from dead.
Postdocs
External
Students