====== 構造最適化 ======

===== Born-Oppenheimer近似 =====

電子系のSchrödinger方程式、
$$
\mathcal{H}_{\mathrm{el}}(\{\vec{R}\})\Psi(\vec{r},\{\vec{R}\})=E(\{\vec{R}\})\Psi(\vec{r},\{\vec{R}\})
$$
のエネルギー固有値$E(\{\vec{R}\})$は原子核の運動に関するポテンシャルとしてはらたく。

したがって、原子核の運動に関するNewtonの運動方程式は
$$
M_{I}\frac{\partial^2}{\partial t^2}\vec{R_I}=-\frac{\partial }{\partial \vec{R}_I}E(\{\vec{R}\})\equiv \vec{F}_{I}
$$

$I$番目の原子核に働く力$\vec{F}_{I}$は
$$
\frac{\partial }{\partial \vec{R}_I}E(\{\vec{R}\})=\left\langle\Psi(\vec{r},\{\vec{R}\})\middle|\frac{\partial }{\partial \vec{R}_I}\mathcal{H}_{\mathrm{el}}(\{\vec{R}\})\middle|\Psi(\vec{r},\{\vec{R}\})\right\rangle
$$
から求めることができる。これをHellmann–Feynman力という。

二階微分のNewtonの運動方程式の代わりに、
$$
\gamma_{I}\frac{\partial}{\partial t}\vec{R_I}= \vec{F}_{I}
$$
のような一階微分の方程式を逐次的に解くことにより、Hellmann–Feynman力がゼロになるような原子配列を得ることができる。
===== CO =====
<code python>
from ase import Atom, Atoms
d = 1.1
atom1 = Atom('C', (0, 0, 0))
atom2 = Atom('O', (d, 0, 0))
molecule1 = Atoms([atom1, atom2], cell = (6, 6, 6), pbc = True)

from ase.calculators.jacapo import Jacapo
solver1 = Jacapo('sample.nc', nbands = 8, ft = 0.01, atoms = molecule1)

from ase.optimize import QuasiNewton
dynamics1 = QuasiNewton(molecule1, trajectory = 'sample.traj')
dynamics1.run(fmax = 0.05)
</code>


練習
  - CO<sub>2</sub>
  - H<sub>2</sub>O
  - NH<sub>3</sub>
  - CH<sub>4</sub>

価電子の数
  - H: 1
  - C: 4
  - N: 5
  - O: 6
  


===== 演習 =====




  - H<sub>2</sub>OのOを固定し二つのHを動かす計算をする<code python>
a1 = (5.0, 0.0, 0.0)
a2 = (0.0, 5.0, 0.0)
a3 = (0.0, 0.0, 5.0)
v1 = (0.50, 0.50, 0.50)
v2 = (0.72, 0.50, 0.50)
v3 = (0.50, 0.72, 0.50)

from ase import Atom, Atoms
p1 = Atom('O', v1, tag = 1) # <--
p2 = Atom('H', v2)
p3 = Atom('H', v3)
molecule1 = Atoms([p1, p2, p3], pbc = True)
molecule1.set_cell([a1, a2, a3], scale_atoms=True)

from ase.constraints import FixAtoms
fixatom1 = [atomx.tag == 1 for atomx in molecule1]
molecule1.set_constraint(FixAtoms(mask = fixatom1))

from ase.calculators.jacapo import Jacapo
solver1 = Jacapo('water.nc', nbands = 8, ft = 0.01, atoms = molecule1)

from ase.optimize import QuasiNewton
dynamics1 = QuasiNewton(molecule1, trajectory = 'water.traj')
dynamics1.run(fmax = 0.04)
print molecule1.get_positions()
print molecule1.get_forces()
</code>  Stay Alive設定をOnにするとどうなるかも確かめよ<code python>
solver1 = Jacapo('water.nc', nbands = 8, stay_alive = True, ft = 0.01, atoms = molecule1)
</code>
  - trajファイルをxyzファイルに変換する<code sh>
$ traj2xyz water.traj water.xyz</code> traj2xyzが無かった場合は[[:script:traj2xyz|traj2xyz]]をつくる。 <code sh>
$ chmod a+x traj2xyz
$ ./traj2xyz water.traj water.xyz</code>
  - water.xyzをiMolかJmolで表示させる。

===== 演習問題 =====
  - 曲がった状態からはじめて、CO<sub>2</sub>の最適化構造を求めよ。
  - 適当な構造からはじめてSiH<sub>4</sub>の最適化構造を求めよ。
  - 上記の問題について、最適化構造に至る過程をアニメーションにせよ。