====== Jacapoチュートリアル ======

===== パラメータの設定 =====
==== 最も簡単な例 ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
atom1 = Atom('C')
atom2 = Atom('O', (1.1, 0, 0))
box1  = Atoms([atom1, atom2], cell = (5, 3, 3), pbc = True)
calc1 = Jacapo(nbands = 6, atoms = box1)
calc1.calculate()
</code>

==== 様々な計算パラメータの設定法 ====
<code python>
from ase import Atom
atom1 = Atom('C', (2, 2, 2))
atom2 = Atom()
atom2.number = 8
(x, y, z) = atom1.position
atom2.position = (x + 1.1, y, z)
from ase import Atoms
box1  = Atoms()
box1.set_cell((5, 3, 3))
box1.set_pbc(True)
box1.append(atom1)
box1.append(atom2)
from ase.calculators.jacapo import Jacapo
ncfile1 = 'sample.nc'
params1 = dict(nbands = 8)
params1.update(stay_alive = False)
params1.update(pw = 200)
params1.update(dw = 400)
params1.update(kpts = (2, 2, 2))
params1.update(ft = 0.01)
params1.update(symmetry = True)
params1.update(atoms = box1)
calc1 = Jacapo(ncfile1, **params1)
calc1.calculate()
print calc1.get_potential_energy()
</code>

===== ncファイルの利用法 =====

==== 使用するncファイルの作成 ====

<code python>
a = 2.0
d = 1.1
from ase import Atom
atom1 = Atom('C', (a, a, a))
atom2 = Atom('O', (a + d, a, a))
from ase import Atoms
box1 = Atoms(cell = (2*a, 2*a, 2*a), pbc = True)
box1.append(atom1)
box1.append(atom2)
from ase.calculators.jacapo import Jacapo
ncfile1 = 'sample.nc'
params1 = dict(nbands = 8)
params1.update(ft = 0.01)
params1.update(symmetry = True)
params1.update(atoms = box1)
calc1 = Jacapo(ncfile1, **params1)
calc1.calculate()
</code>

==== 計算済みncファイルの再利用 ====

<code python>
from ase.calculators.jacapo import read
model1, solver1 = read('sample.nc')
print 'E = ', model1.get_potential_energy(), '[eV]'
print model1.get_forces()
eng = solver1.get_eigenvalues()
occ = solver1.get_occupation_numbers()
for i in range(len(eng)):
  print 'n(%6.2f) = %3.1f'%(eng[i], occ[i])
</code>

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

==== 半手動でCO間距離を変化させる ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
for d in [1.05, 1.10, 1.15, 1.20, 1.25]:
  atom1 = Atom('C', (0, 0, 0))
  atom2 = Atom('O', (d, 0, 0))
  molecule1 = Atoms([atom1, atom2], cell = (6, 6, 6), pbc = True)
  ncfile1 = 'sample_%1.2f.nc' % d
  solver1 = Jacapo(ncfile1, nbands = 8, ft = 0.01, atoms = molecule1)
  print 'd = %1.2f ang' % d
  print 'energy = %f eV' % molecule1.get_potential_energy()
  print 'forces = (eV/ang)\n', molecule1.get_forces()
  print '-----'
  solver1.strip()
</code>

==== グラフ用ファイルの作成 ====

<code python>
from ase.calculators.jacapo import Jacapo
bondlenths =  [1.05, 1.10, 1.15, 1.20, 1.25]
ncfiles = ['sample_%1.2f.nc' % d for d in bondlenths]
solvers = [Jacapo(f) for f in ncfiles]
energies = [s.get_potential_energy() for s in solvers]
print "bondlenth energy"
emin = min(energies)
for i in range(len(bondlenths)):
  print bondlenths[i], energies[i] - emin
</code>
==== 全自動で構造最適化（QuasiNewton版） ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
from ase.optimize import QuasiNewton
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)
solver1 = Jacapo('sample.nc', nbands = 8, ft = 0.01, atoms = molecule1)
dynamics1 = QuasiNewton(molecule1, trajectory = 'sample.traj')
dynamics1.run(fmax = 0.05)

</code>

==== 全自動で構造最適化（BFGS版） ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
from ase.optimize import BFGS
d = 1.15
atom1 = Atom('C', (0, 0, 0))
atom2 = Atom('O', (d, 0, 0))
molecule1 = Atoms([atom1, atom2], cell = (6, 6, 6), pbc = True)
solver1 = Jacapo('sample.nc', nbands = 8, ft = 0.01, atoms = molecule1)
dynamics1 = BFGS(molecule1, trajectory = 'sample.traj')
dynamics1.run(fmax = 0.04)
print molecule1.get_positions()
print molecule1.get_forces()
</code>

==== 拘束条件付き構造最適化 ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
d = 1.20
atom1 = Atom('C', (0, 0, 0), tag = 1)
atom2 = Atom('O', (d, 0, 0))
molecule1 = Atoms([atom1, atom2], cell = (6, 6, 6), pbc = True)
fixatom1 = [atomx.tag == 1 for atomx in molecule1]
molecule1.set_constraint(FixAtoms(mask = fixatom1))
solver1 = Jacapo('sample.nc', nbands = 8, stay_alive = True, ft = 0.01, atoms = molecule1)
dynamics1 = QuasiNewton(molecule1, trajectory = 'sample.traj')
dynamics1.run(fmax = 0.04)
print molecule1.get_positions()
print molecule1.get_forces()
</code>

==== 拘束条件の設定法 ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
from ase.constraints import FixAtoms
from ase.optimize import BFGS

d = 1.20
atom1 = Atom('C', (0, 0, 0), tag = 1)
atom2 = Atom('O', (d, 0, 0))
molecule1 = Atoms([atom1, atom2], cell = (6, 6, 6), pbc = True)
#fix1 = FixAtoms(mask = [atomx.symbol == 'C' for atomx in molecule1])
#fix1 = FixAtoms(indices = [atomx.index for atomx in molecule1 if atomx.symbol == 'C'])
#fix1 = FixAtoms(indices = [atomx.index for atomx in molecule1 if atomx.tag == 1])
fix1 = FixAtoms(mask = [atomx.tag == 1 for atomx in molecule1])
molecule1.set_constraint(fix1)

solver1 = Jacapo('sample.nc', nbands = 8, ft = 0.01, atoms = molecule1)
dynamics1 = BFGS(molecule1, trajectory = 'sample.traj')
dynamics1.run(fmax = 0.04)
print molecule1.get_positions()
print molecule1.get_forces()
</code>

==== trajファイルをxyzファイに変換 ====

<code python>
$ traj2xyz sample.traj sample.xyz
</code>

xyzファイルはiMol/jmolで可視化（アニメーション化）できる。

===== 磁気モーメント（スピン分極）のあるときの計算 =====
==== 強磁性鉄の計算 ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
a = 2.87
atom1 = Atom('Fe', (0.0, 0.0, 0.0), magmom = 1.0)
atom2 = Atom('Fe', (0.5, 0.5, 0.5), magmom = 1.0) 
crystal1 = Atoms([atom1, atom2], pbc = True)
crystal1.set_cell((a, a, a), scale_atoms = True)
para1 = dict(spinpol = True)
para1.update(kpts=(8, 8, 8))
para1.update(pw = 300)
para1.update(dw = 500)
para1.update(symmetry = True)
para1.update(atoms = crystal1)
solver1 = Jacapo('ferro.nc', nbands = 16, **para1)
print solver1.get_potential_energy()
</code>

==== 非磁性鉄の計算 ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
a = 2.87
atom1 = Atom('Fe', (0.0, 0.0, 0.0), magmom = 0.0)
atom2 = Atom('Fe', (0.5, 0.5, 0.5), magmom = 0.0) 
crystal1 = Atoms([atom1, atom2], pbc = True)
crystal1.set_cell((a, a, a), scale_atoms = True)
para1 = dict(spinpol = False)
para1.update(kpts=(8, 8, 8))
para1.update(pw = 300)
para1.update(dw = 500)
para1.update(symmetry = True)
para1.update(atoms = crystal1)
solver1 = Jacapo('nonmag.nc', nbands = 16, **para1)
print solver1.get_potential_energy()
</code>

==== 反強磁性鉄の計算 ====

<code python>
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
a = 2.87
atom1 = Atom('Fe', (0.0, 0.0, 0.0), magmom = 1.0)
atom2 = Atom('Fe', (0.5, 0.5, 0.5), magmom = -1.0) 
crystal1 = Atoms([atom1, atom2], pbc = True)
crystal1.set_cell((a, a, a), scale_atoms = True)
para1 = dict(spinpol = True)
para1.update(kpts=(8, 8, 8))
para1.update(pw = 300)
para1.update(dw = 500)
para1.update(symmetry = True)
para1.update(atoms = crystal1)
solver1 = Jacapo('antiferro.nc', nbands = 16, **para1)
print solver1.get_potential_energy()

</code>