Transcript Slides - Indico

Ljupčo Pejov
[email protected]
Institute of Chemistry, Faculty of Natural
Sciences and Mathematics
Ss. Cyril and Methodius University
Skopje, Macedonia
FP7-INFRASTRUCTURES-2010-2
HP-SEE
High-Performance Computing
Infrastructure for South East Europe’s
Research Communities
THE UNDERLYING PHYSICAL LAWS NECESSARY FOR THE
MATHEMATICAL THEORY OF A LARGE PART OF PHYSICS AND THE
WHOLE OF CHEMISTRY ARE THUS COMPLETELY KNOWN, AND THE
DIFFICULTY IS ONLY THAT THE EXACT APPLICATION OF THESE
LAWS LEADS TO EQUATIONS MUCH TOO COMPLICATED
TO BE SOLUBLE.
THE METHODS OF THEORETICAL PHYSICS SHOULD BE APPLICABLE TO
ALL THOSE BRANCHES OF THOUGHT IN WHICH THE ESSENTIAL
FEATURES ARE EXPRESSIBLE WITH NUMBERS.
AS LONG AS QUANTUM MECHANICS IS CORRECT, ALL ISSUES IN PHYSICS
AND CHEMISTRY ARE ACTUALLY PROBLEMS IN APPLIED MATHEMATICS.
THE LEVEL OF EXACTNESS OF A GIVEN SCIENCE IS DIRECTLY
PROPORTIONAL TO THE LEVEL AT WHICH THE PROBLEMS IN THAT
SCIENCE CAN BE SUBJECTED TO COMPUTATION.
GOD USED BEAUTIFUL MATHEMATICS IN CREATING THE WORLD.
-- P. A. M. DIRAC
 Why do we need computational resources in chemistry






(and chemical engineering).
Which e-infrastructures are more suitable for a particular
problem: grid vs. high-performance clusters.
Choosing and defining the problem.
Developing a computational algorithm and strategy.
Writing the codes.
Implementing the computational approach on a given
platform (linking, compiling, using parallel computing
environment…).
Automating the procedure (???).
 HC-MD-QM-CS - Hybrid Classical or Quantum Molecular
Dynamics – Quantum Mechanical Computer Simulation of
Condensed Phases.
 Application description and objectives
 To study the properties of condensed phases, liquids (such as e.g.
solutions of ions and various molecular systems in molecular liquids),
solids (including small molecular systems adsorbed on surfaces),
computer simulations, which will use parallel numerical algorithms
will be carried out.
 The overall objective of the work is to develop a novel general method
for computation of complex in-liquid properties of the system, with
potential applicability in the field of biomedical sciences, catalysis,
etc.
 Scientific impact:
 These studies are of high fundamental significance, concerning the
properties of condensed phases and influence thereof on various
molecular species
 HPC implementation:

All parts of this application require vast computational efforts, in
particular the quantum molecular dynamics simulations. Though
computations required for certain sequences of the complex hybrid
algorithms can be carried out on standard PC clusters, the more complex
and demanding ones require low-latency parallel environment.
 Achievements:
 Efficient implementation of a genetic algorithm to optimize the interaction





potential energy parameters of liquid CCl4
Efficient implementation of a complex hybrid quantum mechanical (ab initio)
molecular dynamics – quantum mechanical methodology used for simulation of
complex condensed phase systems. The methodology is inherently a multistep one,
and certain steps are excellently scalable on parallel high-performance clusters.
Implementation of hybrid QMD-TDDFT methodology to explain the origin of the
solid-state thermochromism and thermal fatigue of polycyclic overcrowded enes.
Implementation of a hybrid CPMD-QMel-QMnuc methodology to compute the
vibrational spectrum of aqueous hydroxide anion.
Application of molecular dynamics simulations to confirm the experimental
evidence for the existence of dangling OH bonds in hydrophobic hydration shells.
Application of sequential MD-QMel-QMnuc methodology to compute the
vibrational spectra of the first-shell water molecules in ionic water solutions.
Vibrations of aqueous hydroxide anion
Summary of experimental (Raman) and theoretical (power spectrum or
DOS) OH and H2O frequencies from the literature and from our work. H2Osol means water peak in the solution.
 Method
n(OH-) - the anharmonic OH stretching vibrational frequency, i.e. the 0 → 1
vibrational transition
Δn (OH-) - the gas-to-solution frequency shifts, i.e.
Δn (OH-) = n (OH- in solution) - n (isolated OH- ion)
computed from the positions of the peak maxima.
We use a sequential Car-Parrinello Molecular Dynamics (CPMD) +Quantum
chemical (QCelec) +Quantum mechanical (QMnucl) computational approach.
The MD box used in the current study consists of 44 H2O + 2 Na+ + 2 OH(corresponding to a 2.5 molal concentration). Our calculated spectrum will be
compared to the experimental Raman study of a 2.0 molal NaOH(aq) solution,
which has a reported peak maximum at 3625 cm-1 in the region of interest.
We will use the label O* for the hydroxide ion's oxygen atom, and H* for its
hydrogen atom, contrary to Ref.
 Step 1, aqueous solution simulations.
Car-Parinello Molecular Dynamics (CPMD) simulations
of aqueous NaOH solutions.
Hˆ  
i j
ˆ
Pi
2
2M i

1
40

i j
 
 
Zi Z j e2
   U elec Ri ; r
Ri  Ri


d Ri
Mi
  Ri U Ri t 
2
dt

2
d Ri
ˆ 
 min  H
Mi


0
e
0
2
Ri
dt
2




LCP



 N  N
 2
N
1
KS
 i    M I RI    
i 
 i  E  i , R
R , R ,  i , 
I 2
i

E KS

M I RI       ij   i  j
RI
RI
ij
KS

E
 i  

   ij  j
 i
j

 Leap-frog algorithm


1
1


vi (t  t )  vi (t  t )
dvi (t )
Fi
2
2


dt
t
mi

Fi


1
1
vi (t  t )  vi (t  t ) 
t
2
2
mi
1
t  t
2

v
i,
1
2

v
i,
1
2

Fi

t
mi
t
t
1
t
2



1
1
vi (t  t )  vi (t  t )  ai t
2
2



ri (t  t )  vi t  ri (t )
 Verlet algorithm
1
1
2
r t  t   r (t )  v(t )t  a(t )t  b(t )t 3  O(t 4 )
2
6
1
1
2
r t  t   r (t )  v(t )t  a(t )t  b(t )t 3  O(t 4 )
2
6
r t  t   2r (t )  r (t  t )  a(t )t 2  O(t 4 )
Density functional theory

 
 
*
 r     r  ri ψd1d 2  d N    r  ri  
 

 r1   N  1 ,  2 ,...,  N  N
2
d 1d 2 d N
 
  r dr  N
Hˆ   E
a
ˆ
E   H N
a
N
Hˆ 0  Hˆ e  Uˆ ej  Uˆ ee
2

Hˆ e  
2 m0
Uˆ ee
2
ke

2
N
N
2

 el
l 1
N

l 1 k 1
N
1
 
rl  rk
n
Zs
2
ˆ
U ej  ke   
l 1 s 1 rl  Rs

U ej   uˆ (r )
N
l 1
n
Zs

2
uˆ (r )  ke   
s 1 rl  R s
2

E   Na 
2m 0
2
ke
2


el 
2
l 1
N
N
N

l 1 k 1
2
  

E   vˆ ext r  r dr   Na 
2m 0
N
 a
1
    uˆ (r ) N
rl  rk
l 1
2
e
2


ei 
4 0
i 1
N
N
1 a
N

i  j rij
  
E   T    U ej    U ee     vˆr  r dr  FHK  
FHK    T    U ee  
  
E    vˆext r  r dr  Ts  Vclass    Exc  
 2

  vˆext  vˆclass  vˆ xc  i   i i

 2m0

Exc    Ex    Ec  
1
3
1
3 3  3 
LDA
E x        dr
4  
EcVW N
2



bx
x

x
2b  2 x0 
A
x
2b
Q
Q 
1
1
0 
0
 ln

tan

ln

tan



2  X x  Q
2 x  b X x0  
X x 
Q
2 x  b 
 4 
rs  

 3 
1
2
s
xr

1
3
X ( x)  x  bx  c
2

Q  4c  b
E
GGA
xc

1
2 2


2
    dr  xc  ,  ,  

 Each CPMD simulation consists of two phases:
wavefunction optimization and
molecular dynamics simulation.
As the optimization runs involve an iterative procedure that needs to
converge, the number of iterations required to achieve final
convergence being strongly dependent on the particular architecture,
compilers and compilation parameters etc., this phase is strongly
platform – dependent and not so suitable for benchmarking.
The molecular dynamics phase. To demonstrate the scalability of the
approach, in Fig. 1 variation of the wall-clock computational time
required to carry out an MD simulation of a modest-size water cluster
(consisting of 32 water molecules) is plotted against the number of
computing processes (processors/cores). Fig. 2 shows the same data,
where both axes are logarithmic (log-log plot). The parallelization has
been achieved by the MPI paradigm.
8000
7000
6000
t/s
5000
4000
3000
2000
1000
0
0
5
10
15
20
n
Fig. 1.
25
30
35
t/s
10000
1000
100
1
10
n
Fig. 2.
100
 We selected 76 snapshots from the CPMD simulation,
taken at regular intervals from the 10 ps long
production run, i.e. 152 OH- candidates for a
vibrational analysis.
 We then excluded from the analysis all configurations
where proton transfers (or good transfer attempts)
were found to occur. We classified a snapshot as a
proton transfer case along the Ow - - -O* coordinate, if
for
, at least one of the following three
criteria is fulfilled: (i) r1 > 1.2 Å, (ii) r2 < 1.4 Å, (iii)
|δ| = |r1 − r2 | < 0.1 Å.
 The remaining number of cases to be analysed was
130, which we thus believe to represent geometries
sufficiently far from a proton transfer.
 Step 2, construction of system to be used in the
Quantum Chemical (QC) calculations. From each of the
collected snapshots, clusters were extracted, for which we
performed QC calculations in Step 3. A dividing plane
between the O* and H* sides of the ion is defined in Fig. 3.
Each "QM cluster" was composed of a central OH- ion, all
water molecules residing within an R(Ow • • • O*) distance
of less than 3.5 Å on the O* side, and all water molecules
residing within an R(Ow • • • H*) distance of less than 4.5 Å
on the H* side. These cut-offs were based on the
information from the radial pair distribution functions
(rdfs).
Definition of regions I, II and III constructed with
the help of a dividing plane through the center of
the O*–H* bond of the hydroxide ion. Region I
contains all water molecules within a R(O*–Ow)
distance of 3.5 A. Regions II and III contain all
water molecules within a R(H*–Ow) distance of 4.5
A.
 Step 3, QC calculations of the potential energy
surface (PES). One-dimensional potential energy
curves ΔE(r(OH-)) (or U (r(OH-)) were calculated in
the interval 0.7 Å < r(OH-) < 1.5 Å with an increment
of 0.015 Å. The B3LYP/6-31G++(d,p) method with
basis-set superposition error (BSSE) correction
according to the Counterpoise method.
p12
p 22
1
2
H ( x1 , x2 , p1 , p 2 ) 

 k x2  x1  x0 
2m1 2m2 2
r  x2  x1  x0
m1 x1  m2 x2
s
m1  m2
H
  p i
qi
p1  m1 x1

r 
p2  m2 x 2

1
2
2


T  m1 x1  m2 x 2
2
1  m2 x 2
m
x
1
x 2  x1 s 
m1  m2
m2
x1  s 
r
m1  m2 
x 2  s 
m1
r
m1  m2 
 Step 4, vibrational Quantum Mechanical (QM)
calculations. The vibrational energy levels were
calculated quantum-mechanically from the onedimensional potential energy curves using the discrete
variable representation (DVR)
       i     i     i    i
i
i
     k i  k  k i    k i  k 
k
i
k 
2
N
i
  F
k

k
i
0 
0 
 

 
 n  1  n - th row

 
0 
0 
 

 
H in  F TF  V  n
1
i
 Step 5, analysis of frequency vs. field correlations. We
will discuss the frequency shifts in relation to the electric
field that the full hydration shell, or selected parts of it,
generate over the molecule. Here we will probe the electric
field at the equilibrium H* position for each snapshot (the
equilibrium position is known from Step 3) and, moreover,
only probe the component along the O*-H* bond. We
denote this F// @ H*. The // subscript means "along the
vibrational coordinate", i.e. here along the OH bond. We
define the electric field as positive along the molecular axis
if it is oriented as if there were positive charges on the
oxygen side of the molecule and negative charges on the H
side
+
–
+
+
[ O-H]–
–
–
positive field direction
Graphical representation of the positions and orientations of the water
molecules relative to the ion in the CPMD simulation. Lower panel: each dot
represents the position of the center-of-mass of a molecule. Upper panels:
Radial distribution functions for the three regions I, II, and III. The curves have
been scaled according to the respective volumes of regions I, II and III to reflect
the number of water molecules in each region, i.e. they take the solid angle of
the region into account.
Sample potential energy curve for E(rOH) vs. rOH for one of the OH ions
in one of the snapshots from the CPMD simulation and potential
energy curve for the isolated OH ion. The curves have been made to
coincide at their respective minima. The difference between the two
curves is the curve marked Eext, where ext stands for external, i.e. the
ion’s interaction with its surroundings. BSSEcorrectedB3LYP/631G++(d,p) calculations.
Anharmonic OH frequency histograms. (a)
finite QC cluster, (b) periodic QC
calculations.
Conclusions
Which is the best way to model molecules in liquids ?
”Better - in principle”
MD + QM
?
Quantum-dynamical Car-Parrinello
Many-body potentials from ab
initio
QC for some fixed geometry?




ISIS and KALKYL HPC clusters at UPPMAX – Uppsala, Sweden
SEE-grid cluster
The new HPC cluster at FINKI, Skopje, Macedonia
The new HPC cluster at Institute of Chemistry, Skopje,
Macedonia (The FEYNMAN cluster)
THANK YOU FOR YOUR ATTENTION