Transcript Document

South China University of Technology
Optimization for
first principle calculations
Xiaobao Yang
Department of Physics
www.compphys.cn
www.compphys.cn
Review of first principle calculations
ab inito (quantum chemistry)
Density Functional Theory
Exc: LDA GGA(PBE)
Pseudo potential
www.compphys.cn
Lower, more stable
Man fights upwards while water flows downwards.
www.compphys.cn
Inverse Problem and Curse of Dimensionality
www.compphys.cn
Dense packing
S.Torquato
Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. I. Algorithmic details (2005)
www.compphys.cn
www.compphys.cn
Self-assembly of uniform polyhedral silver nanocrystals into
densest packings and exotic superlattices
www.compphys.cn
Root of Equation
www.compphys.cn
www.compphys.cn
www.compphys.cn
Quasi-Newton method
www.compphys.cn
Quasi-Newton method in VASP
www.compphys.cn
Quasi-Newton method in VASP
www.compphys.cn
Steepest descent method
www.compphys.cn
Searching direction
Optimized step
New start point
www.compphys.cn
Searching direction
Optimized step
New start point
New start point
Stop searching
www.compphys.cn
设有二次函数：f(x) = 1/2 (x - x*)TA(x - x*) , 其中A是n×n对称正定矩阵，x*是
一个定点，函数f(x)的等值面 1/2 (x - x*)TA(x - x*) = c 是以x*为中心的椭球面.
由于 ▽f(x*) = A(x - x*) = 0，A正定，因此x*是f(x)的极小点。设x(1)是在某个等
值面上的一点，该等值面在点x(1)处的法向量▽f(x(1)) = A(x(1) - x*)。
设d(1)是这个等值面在x(1)处的一个切向量, d(1)与▽f(x(1))正交，即
d(1)T▽f(x(1)) = 0; 另外, d(2) = x* - x(1), 有 d(1)TAd(2) = 0，
等值面上一点处的切向量与由这一点指向极小点的向量关于A共轭
www.compphys.cn
Conjugate Gradient method
Steepest descent method
Conjugate Gradient
www.compphys.cn
www.compphys.cn
Structural relaxation
To obtain the ground state relaxed geometry of the system.
the equilibrium lattice constants
a given ionic configuration
the forces obtained
these forces are greater
than some minimum tolerance
the ions are moved in the direction of the forces
www.compphys.cn
IBRION-tag in VASP
0 Standard ab-initio molecular dynamics
A Verlet algorithm is used to integrate Newton's equations of
motion. POTIM supplies the timestep in femto seconds.
1 A quasi-Newton algorithm is used to relax the ions into
their instantaneous ground state.
2 A conjugate-gradient algorithm
i)In the first step ions :the direction of the steepest descent. The
energy and the forces are recalculated. ii) Interpolation of the
change of the total energy and of the forces, then a corrector
step. iii) After the corrector step the forces and energy are
recalculated and it is checked whether the forces contain a
significant component parallel to the previous search direction.
1 initial position 2 trial step 3 corrector step, i.e. positions
corresponding to anticipated minimum 4 trial step
www.compphys.cn
5 corrector step …
Relaxation with VASP
www.compphys.cn
Energies vs. Forces
H2 molecular
C6H6 molecular
C
C
C
C
C
C
H
H
H
H
H
H
-2.629
-1.415
-1.493
-2.618
-3.851
-3.904
-2.682
-0.585
-0.544
-2.613
-4.764
-4.683
3.681
2.939
1.580
0.818
1.541
2.981
4.664
3.486
1.065
-0.192
1.084
3.483
0.018
0.026
0.024
-0.011
0.016
-0.033
0.020
-0.047
-0.022
-0.045
-0.040
0.032
www.compphys.cn
heuristic algorithm for global optimization
一个基于直观或经验构造的算法，在可接受的花费
（指计算时间和空间）下给出待解决组合优化问题
每一个实例的一个可行解，该可行解与最优解的
偏离程度不一定事先可以预计.
Simulated Annealing
Genetic Algorithm
Particle Swarm Optimization
www.compphys.cn
No Free Lunch
最优化理论的发展之一是wolpert和
Macerday提出了没有免费的午餐定理
(No Free Lunch，简称NFL)。该定理的
结论是，由于对所有可能函数的相互补偿
，最优化算法的性能是等价的。该定理暗
指，没有其它任何算法能够比搜索空间的
线性列举或者纯随机搜索算法更优。该定
理只是定义在有限的搜索空间，对无限搜
索空间结论是否成立尚不清楚。
www.compphys.cn
Simulated Annealing
s ← s0; e ← E(s)
// Initial state, energy.
sbest ← s; ebest ← e
// Initial "best" solution
k←0
// Energy evaluation count.
while k < kmax and e > emax
// While time left & not good enough:
T ← temperature(k/kmax)
// Temperature calculation.
snew ← neighbour(s)
// Pick some neighbour.
enew ← E(snew)
// Compute its energy.
if P(e, enew, T) > random() then
// Should we move to it?
s ← snew; e ← enew
// Yes, change state.
if e < ebest then
// Is this a new best?
sbest ← snew; ebest ← enew
// Save 'new neighbour' to 'best found'.
k←k+1
// One more evaluation done
return sbest
// Return the best solution found.
www.compphys.cn
Structure of an MD code
www.compphys.cn
Key issues in MD
►Classical Molecular Dynamics
- Potentials to describe the particles in the system
formulism + parameterization
- Algorithm to solve Newton EQ:
accuracy + efficiency
- Analysis of physical properties
►Ab initio Molecular Dynamics
- Coupling the motion of electrons and ions.
www.compphys.cn
Why classical MD works for Atoms?
►Energy View:
Typical Kinetic energy < 0.1 eV, while it is about >1eV to
remove/excite an elec.
►Debroglie wavelength view:
Typical distance between atoms > 1 Å, while the DeBroglie
wavelength is ~ 10-7 Å.
►MD vs. MC:
MC: be convenient for studying the equilibrium properties.
MD: reflects the real process of a system from one microstate to another.
www.compphys.cn
EoM: Verlet method
Equation of Motion of particle i:
Taylor expansion of positions with time
d 2 ri
 a(ri , t )
2
dt
Verlet method
What is the difference with Euler or Euler-Cromer method?
www.compphys.cn
Properties of a dilute gas
► Lennard-Jones potential
V ( r )  4 [(
f k ,i

r
)12  (

r
)6 ]
V
 12
6

 24 [( 13 )  ( 7 )]
rk ,i
r
r
dvi , x
dt
dvi , y
dt
 ai , x ,
 ai , x ,
dxi
 vi , x
dt
dyi
 vi , y
dt
1
ai , x   f k ,i cos  k ,i
m k i
1
ai , y   f k ,i sin  k ,i
m k j
www.compphys.cn
Morse potential
www.compphys.cn
Many body potential
www.compphys.cn
Properties of a dilute gas
► Periodical boundary condition
PBC affects the force
calculation and position update.
► thermodynamic ensembles
Temperature:
3
NkT  K  t 
2
time
1 N
  mvi2  t 
2 i 1
time
Instantaneous temperature (T*):
N
3
1
NkT *  t   K  t    mvi2  t 
2
2 i 1
www.compphys.cn
Structure of an MD code
www.compphys.cn
MD.m
[ pos, vel, acc, seed ] = initialize ( np, nd, box, seed );
[ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );
for step = 1 : step_num
[ force, potential, kinetic ] = compute ( np, nd, pos, vel, mass );
[ pos, vel, acc ] = update ( np, nd, pos, vel, force, acc, mass, dt );
end
Initialize.m
Compute.m
Update.m
www.compphys.cn
Initialize.m

pos(1:nd,1:np) = rand ( nd, np );








for i = 1 : nd
pos(i,1:np) = box(i) * pos(i,1:np);
end
vel(1:nd,1:np) = 0.0;
acc(1:nd,1:np) = 0.0;
return
end
www.compphys.cn
Compute.m






%v(x) = ( sin ( min ( x, PI2 ) ) )**2
%dv(x) = 2.0 * sin ( min ( x, PI2 ) ) * cos
( min ( x, PI2 ) )= sin ( 2.0 * min ( x, PI2 ) )
pot = pot + 0.5 * sum ( sin ( D2 ).^2 );
% accumulate pot. energy
f( :, i) = Ri * ( sin ( 2*D2 ) ./ D )';
% force on particle 'i'
% Compute kinetic energy.
kin = 0.5 * mass * sum ( diag ( vel' * vel ) );
www.compphys.cn
Update.m







%
%
%
x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt * dt
v(t+dt) = v(t) + 0.5 * ( a(t) + a(t+dt) ) * dt
a(t+dt) = f(t) / m
pos(1:nd,1:np) = pos(1:nd,1:np) + vel(1:nd,1:np) * dt ...
+ 0.5 * acc(1:nd,1:np) * dt * dt;
vel(1:nd,1:np) = vel(1:nd,1:np) ...
+ 0.5 * dt * ( f(1:nd,1:np) * rmass + acc(1:nd,1:np) );
acc(1:nd,1:np) = f(1:nd,1:np) * rmass;
www.compphys.cn
Properties of a dilute gas
Trajectories
Speed distribution
► 20 particles in a 10x10 box with PBC
kT 
m
(v x 2  v y 2 ) 
2
www.compphys.cn
The Melting Transition
www.compphys.cn
The Melting Transition
www.compphys.cn
MD: isothermal molecular dynamics
How can we modify the EoM so that they lead to
constant temperature?
Nose-Hoover thermostat
Berendsen’s thermostat
Direct feedback
www.compphys.cn
Genetic Algorithm

适应度(fitness)是借鉴生物个体对环境的适应程度, 而对所
求解问题中的对象设计的一种表征优劣的测度。

以生物细胞中的染色体(chromosome)代表问题中的个体对
象。遗传算法中染色体一般用字符串表示, 而基因也就是字符串
中的一个个字符,用二进制数串作为染色体编码。

遗传算法是通过在种群(population-由若干个染色体组成的
群体)上实施所称的遗传操作,使其不断更新换代而实现对整个参
数空间的搜索。

遗传算法中有三种关于染色体的运算: 选择-复制、交叉和变异,
这三种运算被称为遗传操作或遗传算子(genetic operator)。
www.compphys.cn
选择-复制
选择-复制(selection
reproduction)操作是模拟生物界优胜劣汰的自
然选择法则的一种染色体运算, 就是从种群中选择适应度较高的染色体进行复制,以
生成下一代种群。选择-复制的通常做法是, 对于一个规模为N的种群S,按每个染色体
xi∈S的选择概率P(xi)所决定的选中机会, 分N次从S中随机选定N个染色体, 并进行复
制。 这里的选择概率P(xi)的计算公式为
P( xi ) 
f ( xi )
N
 f (x )
j 1
j
按照这种选择概率定义, 适应度越高的染色体被随机选定的概率就越大,
被选中的次数也就越多, 从而被复制的次数也就越多。
www.compphys.cn
交叉 交叉 (crossover)亦称交换、交配或杂交,就是互
换两个染色体某些位上的基因。例如,设染色体s1=01001011,
s2=10010101, 交换其后4位基因, 即
则得新串s1′=01000101,s2′=10011011。s1′和s2′可以看做是原染
色体s1和s2的子代染色体。
变异 变异(mutation)亦称突变,就是改变染色体某个(些)位
上的基因。例如,把染色体s=11001101的第三位上的0变为1, 则
得到新染色体s′=11101101。
www.compphys.cn
求解区间［0,31］上的二次函数y=x2的最大值


取［0,31］中的整数用5位二进制数作为个体x的基因型编码。
种群规模设定为4,取染色体
s1=01101(13),s2=11000(24),s3=01000(8),
s4=10011(19)组成初始种群S1。
选择-复制
设从区间［0, 1］中产生4个随机数如下:
r1=0.450126, r2=0.110347, r3=0.572496, r4=0.98503
s1’=11000（24）, s2’=01101（13）, s3’=11000（24）, s4’=10011（19）
www.compphys.cn
交叉 设交叉率pc=100%，即S1中的全体染色体都参加交叉运
算。设s1’与s2’配对，s2’与s4’配对。分别交换后两位基因，
得新染色体：
s1’=11000（24）, s2’=01101（13）, s3’=11000（24）, s4’=10011（19）
s1’’=11001（25）, s2’’=01100（12）,
s3’’=11011（27）, s4’’=10000（16）
变异 设变异率pm=0.001。这样，群体S1中共有540.001=0.02
位基因可以变异。0.02位显然不足1位，所以本轮遗传操作不
做变异。
www.compphys.cn
假设这一轮选择-复制操作中，种群S2中的4个染色体都被
选中（因为选择概率毕竟只是一种几率，所以4个染色体恰好
都被选中的情况是存在的），得到群体：
s1’=11001（25）, s2’=01100（12）, s3’=11011（27）, s4’=10000（16）
做交叉运算，让s1’与s2’，s3’与s4’ 分别交换后三位基因，得
s1’’=11100（28）, s2’’=01001（9）, s3’’=11000（24）, s4’’=10011（19）
这一轮仍然不会发生变异。于是，得第三代种群S3：
www.compphys.cn
设这一轮的选择-复制结果为：
s1’=11100（28）, s2’=11100（28）, s3’=11000（24）, s4’=10011（19）
然后，做交叉运算，让s1’与s4’，s2’与s3’ 分别交换后两位基因，得
s1’’=11111（31）, s2’’=11100（28）, s3’’=11000（24）, s4’’=10000（16）
这一轮仍然不会发生变异。于是，得第四代种群S4：
s1=11111（31）, s2=11100（28）, s3=11000（24）, s4=10000（16）
www.compphys.cn
www.compphys.cn
Particle Swarm Optimization
PSO中，每个优化问题的解都是搜索空间中的一只鸟，
我们称之为“粒子”。所有的粒子都有一个由被优化的函
数决定的适应值(fitness value)，每个粒子还有一个速度
决定他们飞翔的方向和距离，然后粒子们就追随当前的最
优粒子在解空间中搜索。
PSO初始化为一群随机粒子(随机解)，然后通过迭代找
到最优解。在每一次迭代中，粒子通过跟踪两个“极值”来
更新自己。一个就是粒子本身所找到的最优解，这个解叫
做个体极值pBest，另一个极值是整个种群目前找到的最
优解，这个极值是全局极值gBest。
www.compphys.cn
www.compphys.cn
www.compphys.cn
粒子更新公式

在找到这两个最优值时，粒子根据如下的公式来更新自己的速度和
新的位置：
其中 ,V 是粒子的速度 , Present 是粒子的当前位置 ,pBest 与
gBest见前面定义。rand ( )是(0 ,1)之间的随机数 ,c1和c2被称
作学习因子。通常 ,c1 = c2 = 2。w 是加权系数（惯性权重） ,
取值在 0. 1到0. 9之间。粒子通过不断学习更新 ,最终飞至解空间
中最优解所在的位置 ,搜索过程结束。最后输出的 gBest 就是全局
最优解。在更新过程中 ,粒子每一维的最大速率限被限制为 Vmax
,如果某一维更新后的速度超过设定的Vmax，那么这一维的速度就
被限定为Vmax。
www.compphys.cn