====== STM像のシミュレーション ======
===== Al(111) with Octave =====
<code python>
# Save this file as Al111.py.
from ase.lattice.surface import fcc111
atoms = fcc111('Al', size=(1, 1, 2))
atoms.center(vacuum=4.0, axis=2)
atoms.pbc=True
from ase.calculators.jacapo import Jacapo
calc = Jacapo(nbands = 10, kpts = (4, 4, 1), atoms = atoms)
energy = atoms.get_potential_energy() 
print energy

from ase.dft.stm import STM
stm = STM(atoms)
z = 8.0
bias = 1.0
c = stm.get_averaged_current(bias, z)
x, y, h = stm.scan(bias, c, repeat=(3, 5))

print '# Created for Octave '
print '# name: x'
print '# type: matrix'
print '# rows: %d'%(len(x))
print '# columns: %d'%(len(x[0]))
for i in range(len(x)):
  for j in range(len(x[0])):
    print x[i][j],
  print
print '# name: y'
print '# type: matrix'
print '# rows: %d'%(len(y))
print '# columns: %d'%(len(y[0]))
for i in range(len(y)):
  for j in range(len(y[0])):
    print y[i][j],
  print
print '# name: z'
print '# type: matrix'
print '# rows: %d'%(len(h))
print '# columns: %d'%(len(h[0]))
for i in range(len(h)):
  for j in range(len(h[0])):
    print h[i][j],
  print

</code>

<code sh>
$ python Al111.py > Al111.data
$ octave
>>> load Al111.data
>>> mesh(x, y, z)
>>> contour(x, y, z)
>>>
</code>


===== from ncファイル to octave =====
<code python>
from ase.calculators.jacapo import read
atoms, calc = read('out.nc')
energy = atoms.get_potential_energy() 
print energy

from ase.dft.stm import STM
stm = STM(atoms)
z = 8.0
bias = 1.0
c = stm.get_averaged_current(bias, z)
x, y, h = stm.scan(bias, c, repeat=(3, 5))

print '# Created for Octave '
print '# name: x'
print '# type: matrix'
print '# rows: %d'%(len(x))
print '# columns: %d'%(len(x[0]))
for i in range(len(x)):
  for j in range(len(x[0])):
    print x[i][j],
  print
print '# name: y'
print '# type: matrix'
print '# rows: %d'%(len(y))
print '# columns: %d'%(len(y[0]))
for i in range(len(y)):
  for j in range(len(y[0])):
    print y[i][j],
  print
print '# name: z'
print '# type: matrix'
print '# rows: %d'%(len(h))
print '# columns: %d'%(len(h[0]))
for i in range(len(h)):
  for j in range(len(h[0])):
    print h[i][j],
  print


</code>


===== Matplotlibによる描画 =====



<code python>
from ase.calculators.jacapo import read
atoms, calc = read('out.nc')

from ase.dft.stm import STM
stm = STM(atoms)
z = 8.0
bias = 1.0
c = stm.get_averaged_current(bias, z)
x, y, h = stm.scan(bias, c, repeat=(3, 5))

import matplotlib.pyplot as plt
plt.gca(aspect='equal')
plt.contourf(x, y, h, 40)
plt.hot()
plt.colorbar()
plt.savefig('2d.png')
plt.show()

</code>

===== コマンドラインからncファイルを指定 =====
<code python>
import sys
argvs = sys.argv
argc = len(argvs)
if (argc < 2):
    print 'Usage: # python %s filename' % argvs[0]
    quit()

from ase.calculators.jacapo import read
atoms, calc = read(argvs[1])
energy = atoms.get_potential_energy() 

from ase.dft.stm import STM
stm = STM(atoms)
z = 8.5
bias = 5.0
c = stm.get_averaged_current(bias, z)
x, y, h = stm.scan(bias, c, repeat=(3, 5))

import matplotlib.pyplot as plt
plt.gca(aspect='equal')
plt.contourf(x, y, h, 40)
plt.hot()
plt.colorbar()
plt.show()

</code>

注意: Matplotlibはsstg3にしかインストールされていない。

[[seminar:jacapo_sample|目次へもどる]]