First-principle molecular dynamics with ultrasoft pseudopotentials: theory, parallel implementation, applications P. Giannozzi Scuola Normale Superiore di Pisa and Democritos National Simulation Center, Italy • Ultrasoft Vanderbilt pseudopotentials in first-principle molecular dynamics • Practical implementation and parallelization • Applications: biologically relevant molecules Simplest model for the myoglobin active site: iron-porphyrin-imidazole complex. Yellow: Fe. Dark gray: C. Blue: N. Light gray: H. Work done in collaboration with F. de Angelis (Perugia) and R. Car (Princeton)
Goal : Modellization of biological molecules (and of their reactions) containing transition metal centers. Extended model of the myoglobin active site, including the 13 surrounding residues in a 8A sphere centered around the iron atom, terminated with NH 2 groups (332 atoms, 902 electrons). Red: O, other atoms as in previous figure. Interest : respiration, photosynthesis, enzymatic catalysis...
Influence of the protein environment Spin density distribution of the quintet state ( S = 2) for the extended model of the myoglobin active site, showing a sizable contribution not localized on the Iron atom
Problems to be addressed: • high-level theoretical description needed, due to the presence of transition metal atoms • complex energy landscape: many possible spin states and local minima, not easy to guess from static (single-point) calculations • large systems (hundreds of atoms) needed in order to accurately describe the effect of surrounding environment A powerful tool: First-Principle Molecular Dynamics • combines forces on atoms from electronic structure with classical Molecular Dynamics simulation machinery • forces on atoms are usually obtained from Density-Functional Theory • a convenient way to perform First-Principle Molecular Dynamics: Car-Parrinello dynamics on both electronic and nuclear degrees of freedom
Density-Functional Theory The energy is a functional of the charge density n ( r ) : � E DF T = F [ n ( r )] + n ( r ) V ( r ) d r that is minimized by the ground-state charge density n ( r ) . Typically solved via Kohn-Sham equations: � � − � 2 2 m ∇ 2 + V ( r ) + V H ( r ) + V xc ( r ) ψ i ( r ) = ǫ i ψ i ( r ) where � | ψ i ( r ) | 2 , n ( r ) = � ψ i | ψ j � = δ ij , i the Hartree potential V H ( r ) is: � n ( r ′ ) V H ( r ) = e 2 | r − r ′ | d r ′ , and V xc ( r ) is a suitable exchange-correlation potential , approximated as a function of n ( r ) and of its gradient |∇ n ( r ) | .
Forces in Density-Functional Theory Forces acting on nuclei have the simple Hellmann-Feynman form: � F I = − ∂E ( { R I } ) n ( r ) ∂V d r + ∂E N = − ∂ R I ∂ R I ∂ R I R I = position of I − th nucleus, E ( { R I } ) = E DF T + E N , E N = nuclear electrostatic energy. Structural minimization : find minimum of E ( { R I } ) , calculated for the electronic ground-state charge density corresponding to nuclei in positions { R I } . Molecular Dynamics : introduce Lagrangian � L = 1 M I ∇ 2 R I − E ( { R I } ) 2 I generating equations of motions: R I = − ∂E ( { R I } ) F I = M I ¨ ∂ R I that can be discretized in time using i.e. Verlet algorithm.
Car-Parrinello Molecular Dynamics Introduce fictitious dynamics on the electronic degrees of freedom ψ i : � � � ψ i ( r ) | 2 + 1 | ˙ M I ∇ 2 L = µ R I − E ( { ψ i } , { R I } ) 2 i I ( µ = fictitious electronic mass), subject to orthonormality constraints � ψ i | ψ j � − δ ij = 0 . Generated equations of motion: � ψ i = − δE µ ¨ + Λ ij ψ j δψ i ij R I = − ∂E F I = M I ¨ ∂ R I where Λ ij = Lagrange multiplier enforcing the constraints. Electrons are initially brought close to the ground state at fixed nuclei. The combined electronic and nuclear dynamics keeps electrons close to the instantaneous ground state.
Practical Implementation A supercell geometry and the corresponding Plane-Wave (PW) basis set are typically used: 1 e i ( k + G ) · r √ � r | k + G � = Ω G = reciprocal lattice vectors of the supercell, Ω = unit cell volume, k = Bloch vector. The Bloch vector can be assumed to be k = 0 and dropped for aperiodic systems (molecules). Expansion of orbitals into PW’s: � ψ n ( G ) 1 e i G · r . √ ψ n ( r ) = Ω G A finite PW set is obtained by taking all PW’s up to a given kinetic energy cutoff E c : � 2 2 m | G | 2 ≤ E c . The basic operation to be performed is the Fourier Transform .
Norm-Conserving Pseudopotentials Electron-ionic core interactions are represented by a nonlocal norm-conserving pseudopotential (NCPP): a soft potential, no core states, no orthonormality wiggles. 3S AE 3P AE 3D AE 3S PS 3P PS 3D PS 7 0 1 2 3 4 5 6 8 r (a.u.) NCPP’s can be recast into the form: � ˆ V ≡ V loc ( r ) + | β n � D nm � β m | nm In many systems, NCPP’s allow accurate calculations with moderate-size ( E c ∼ 10 − 20 Ry ) PW basis sets
Limitations of Norm-Conserving Pseudopotentials NCPP’s are still “hard”and require a large plane-wave basis sets ( E c > 70 Ry ) for: • first-row elements, in particular N, O, F • transition metals, in particular the 3d row: Cr, Mn, Fe, Co, Ni, ... Even if just one atom is “hard”, a high cutoff is required. This translates into large CPU and RAM requirements. Ultrasoft (Vanderbilt) pseudopotentials (USPP) are devised to overcome such a problem: 3d pseudo- and all-electron orbitals for Cu (Laasonen et al, Phys. Rev. B 47, 10142 (1993))
Ultrasoft pseudopotentials � ˆ V US ≡ V loc ( r ) + D lm | β l �� β m | lm Charge density with USPP: � � � | ψ i ( r ) | 2 + n ( r ) = � ψ i | β l � Q lm ( r ) � β m | ψ i � i i lm where the Q lm (“augmentation charges”) are: l ( r ) φ m ( r ) − � l ( r ) � Q lm ( r ) = φ ∗ φ ∗ φ m ( r ) | β l � are “projectors” | φ l � are atomic states (not necessarily bound) | � φ l � are pseudo-waves (coinciding with | φ l � beyond some “core radius”) In practical USPP, the Q lm ( r ) are pseudized . Orthonormality with USPP: � � ψ ∗ � ψ i | S | ψ j � = i ( r ) ψ j ( r ) d r + � ψ i | β l � q lm � β m | ψ j � = δ ij lm � where q lm = Q lm ( r ) d r
Ultrasoft pseudopotentials and PAW Projector Augmented Waves (PAW) method: P. E. Bl¨ ochl, PRB 50 , 17953 (1994) A linear transformation � T connects “true” orbitals | ψ i � to “pseudo” orbitals | � ψ i � : � � � | ψ i � = � T | � ψ i � = | � | φ l � − | � � β l | � ψ i � + φ l � ψ i � l where | φ l � = “true” atomic states, | � φ l � = pseudo-waves and � β l | � φ m � = δ lm ⇒ � T | � φ l � = | φ l � . The pseudo-orbitals are the variational parameters of the calculation. Assuming that in the core region: � | � | � φ l �� β l | � ψ i � ≃ ψ i � l we recover the USPP expression for the charge density n ( r ) . The PAW procedure can be used to reconstruct all-electron orbitals from pseudo-orbitals
Additional Problems in a USPP-PW calculation • Augmentation term in the charge density (may require a higher cutoff than the “soft” part) � � n ( r ) = n soft ( r ) + � ψ i | β l � Q lm ( r ) � β m | ψ i � i lm • Overlap Matrix in orthonormality constraints � ψ i = − δE µ ¨ + Λ ij Sψ j δψ i ij • Additional terms in the forces, plus a term coming from orthonormality constraints � F I = − ∂E Λ ij � ψ i | ∂S + | ψ j � ∂ R I ∂ R I ij
Implementation Solutions • Simplified iterative orthonormalization procedure • Double FFT grid • “Box” grid for augmentation charge - exploits the localized character of Q lm ( r ) in real space Requirements for an effective parallelization • Load balancing : all processes should do more or less the same work • Communication : slower than calculation, inter-process communications must be kept to a minimum • Scalability : work done by N processes should go like 1 /N Parallelization issues • Parallel 3-dimensional FFT • “Box” grid parallelization Requirements: Load balancing , minimal inter-process communications
Basic ingredients of a USPP-PW calculation • Scalar products and various linear algebra: � β ∗ m ( G ′ ) ψ i ( G ′ ) | β l �� β m | ψ i � ⇒ β l ( G ) G ′ • Fast Fourier-Transform (FFT): � � f ( G ) e i G · r f ( r ) = G � 1 � f ( r ) e − i G · r d r f ( G ) = Ω Discretization: � N 1 e 2 πi mj f ( lmn ) e 2 πi li N 2 e 2 πi nk � f ( ijk ) = N 3 lmn � 1 N 1 e − 2 πi mj f ( ijk ) e − 2 πi li N 2 e − 2 πi nk � f ( lmn ) = N 3 N 1 N 2 N 3 ijk where = ( l − 1) h 1 + ( m − 1) h 2 + ( n − 1) h 3 G lmn i − 1 a 1 + j − 1 a 2 + k − 1 = r ijk a 3 N 1 N 2 N 3
Parallel 3D FFT G -space basis sets: Orbital basis set: � 2 G 2 • 2 m ≤ E cut Soft-charge basis set: � 2 G 2 • 2 m ≤ 4 E cut Augmented-charge basis set: � 2 G 2 • 2 m ≤ d × E cut , d ≥ 4 . R - and G -space FFT grids: • “Coarse” grid for orbitals and soft charge density • “Dense” grid for augmented charge density Requisites for good load balancing: • Number of G -vectors in all basis sets approx. the same on all processors • Number of FFT operations in both grids approx. the same on all processors
Recommend
More recommend
