電子系の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力がゼロになるような原子配列を得ることができる。

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)

練習

1. CO₂
2. H₂O
3. NH₃
4. CH₄

価電子の数

1. H: 1
2. C: 4
3. N: 5
4. O: 6

1. H₂OのOを固定し二つのHを動かす計算をする

   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()

   Stay Alive設定をOnにするとどうなるかも確かめよ

   solver1 = Jacapo('water.nc', nbands = 8, stay_alive = True, ft = 0.01, atoms = molecule1)

2. trajファイルをxyzファイルに変換する

   $ traj2xyz water.traj water.xyz

   traj2xyzが無かった場合はtraj2xyzをつくる。

   $ chmod a+x traj2xyz
   $ ./traj2xyz water.traj water.xyz

3. water.xyzをiMolかJmolで表示させる。

1. 曲がった状態からはじめて、CO₂の最適化構造を求めよ。
2. 適当な構造からはじめてSiH₄の最適化構造を求めよ。
3. 上記の問題について、最適化構造に至る過程をアニメーションにせよ。