South China University of Technology

1. South China University of Technology Optimization for first principle calculations Xiaobao Yang Department of Physics www.compphys.cn

2. Review of first principle calculations ab inito (quantum chemistry) Density Functional Theory Exc: LDAGGA(PBE) Pseudo potential

3. Lower, more stable Man fights upwards while water flows downwards.

4. Inverse Problem and Curse of Dimensionality

5. Dense packing S.Torquato Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. I. Algorithmic details (2005)

7. Self-assembly of uniform polyhedral silver nanocrystals into densest packings and exotic superlattices

8. Root of Equation

11. Quasi-Newton method

12. Quasi-Newton methodin VASP

13. Quasi-Newton methodin VASP

14. Steepest descent method

15. Searching direction New start point Optimized step

16. Searching direction New start point New start point Optimized step Stop searching

17. 设有二次函数：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共轭

18. Conjugate Gradient method Steepest descent method Conjugate Gradient

20. 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

21. 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 5 corrector step …

22. Relaxation with VASP

23. Energies vs. Forces H2 molecular C -2.629 3.681 0.018 C -1.415 2.939 0.026 C -1.493 1.580 0.024 C -2.618 0.818 -0.011 C -3.851 1.541 0.016 C -3.904 2.981 -0.033 H -2.682 4.664 0.020 H -0.585 3.486 -0.047 H -0.544 1.065 -0.022 H -2.613 -0.192 -0.045 H -4.764 1.084 -0.040 H -4.683 3.483 0.032 C6H6 molecular

24. heuristic algorithm for global optimization 一个基于直观或经验构造的算法，在可接受的花费 （指计算时间和空间）下给出待解决组合优化问题 每一个实例的一个可行解，该可行解与最优解的 偏离程度不一定事先可以预计. Simulated Annealing Genetic Algorithm Particle Swarm Optimization

25. No Free Lunch 最优化理论的发展之一是wolpert和Macerday提出了没有免费的午餐定理(No Free Lunch，简称NFL)。该定理的结论是，由于对所有可能函数的相互补偿，最优化算法的性能是等价的。该定理暗指，没有其它任何算法能够比搜索空间的线性列举或者纯随机搜索算法更优。该定理只是定义在有限的搜索空间，对无限搜索空间结论是否成立尚不清楚。

26. SimulatedAnnealing 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.

27. Structure of an MD code

28. Key issues in MD • Classical Molecular Dynamics - Potentials to describe the particles in the system formulism + parameterization formulism + parameterization - Algorithm to solve Newton EQ: accuracy + efficiency accuracy + efficiency - Analysis of physical properties • Ab initio Molecular Dynamics - Coupling the motion of electrons and ions.

29. 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.

30. EoM: Verlet method Equation of Motion of particle i: Taylor expansion of positions with time Verlet method What is the difference with Euler or Euler-Cromer method?

31. Properties of a dilute gas • Lennard-Jones potential

32. Morse potential

33. Many body potential

34. Properties of a dilute gas • Periodical boundary condition PBC affects the force calculation and position update. • thermodynamic ensembles Temperature: Instantaneous temperature (T*):

35. Structure of an MD code

36. 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

37. 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

38. 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 ) );

39. 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;

40. Properties of a dilute gas Trajectories Speed distribution • 20 particles in a 10x10 box with PBC

41. The Melting Transition

42. The Melting Transition

43. 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

44. Genetic Algorithm • 适应度(fitness)是借鉴生物个体对环境的适应程度, 而对所求解问题中的对象设计的一种表征优劣的测度。 • 以生物细胞中的染色体(chromosome)代表问题中的个体对象。遗传算法中染色体一般用字符串表示, 而基因也就是字符串中的一个个字符,用二进制数串作为染色体编码。 • 遗传算法是通过在种群(population-由若干个染色体组成的群体)上实施所称的遗传操作,使其不断更新换代而实现对整个参数空间的搜索。 • 遗传算法中有三种关于染色体的运算: 选择-复制、交叉和变异,这三种运算被称为遗传操作或遗传算子(genetic operator)。

45. 　　选择-复制　选择-复制(selectionreproduction)操作是模拟生物界优胜劣汰的自然选择法则的一种染色体运算, 就是从种群中选择适应度较高的染色体进行复制,以生成下一代种群。选择-复制的通常做法是, 对于一个规模为N的种群S,按每个染色体xi∈S的选择概率P(xi)所决定的选中机会, 分N次从S中随机选定N个染色体, 并进行复制。 这里的选择概率P(xi)的计算公式为 按照这种选择概率定义, 适应度越高的染色体被随机选定的概率就越大, 被选中的次数也就越多, 从而被复制的次数也就越多。

46. 　　交叉 交叉 (crossover)亦称交换、交配或杂交,就是互换两个染色体某些位上的基因。例如,设染色体s1=01001011, s2=10010101, 交换其后4位基因, 即 则得新串s1′=01000101,s2′=10011011。s1′和s2′可以看做是原染色体s1和s2的子代染色体。 变异 变异(mutation)亦称突变,就是改变染色体某个(些)位上的基因。例如,把染色体s=11001101的第三位上的0变为1, 则得到新染色体s′=11101101。

47. 求解区间［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）

48. 　交叉 设交叉率pc=100%，即S1中的全体染色体都参加交叉运算。设s1’与s2’配对，s2’与s4’配对。分别交换后两位基因，得新染色体： s1’’=11001（25）, s2’’=01100（12）, s3’’=11011（27）, s4’’=10000（16） 　变异 设变异率pm=0.001。这样，群体S1中共有540.001=0.02位基因可以变异。0.02位显然不足1位，所以本轮遗传操作不做变异。 s1’=11000（24）, s2’=01101（13）, s3’=11000（24）, s4’=10011（19）

49. 假设这一轮选择-复制操作中，种群S2中的4个染色体都被选中（因为选择概率毕竟只是一种几率，所以4个染色体恰好都被选中的情况是存在的），得到群体：假设这一轮选择-复制操作中，种群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：

50. 设这一轮的选择-复制结果为： 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）