====== Band計算 ======
===== スクラッチから計算する =====

<code python>
""" Make a bandstructure plot using the Harris functional. 
"""
import sys
from Dacapo import Dacapo
from ASE import Atom,ListOfAtoms
from ASE.Utilities.List import List
import Numeric as num
import os

bulk = ListOfAtoms([Atom('Cu', (0,    0,     0))] )
a0_fcc = 3.61
b = a0_fcc/2.
bulk.SetUnitCell([(0, b, b), (b, 0, b), (b, b, 0)])

calc = Dacapo()
calc.SetBZKPoints((10,10,10))        # set the k-points Monkhorst-Pack 
calc.SetPlaneWaveCutoff(340)         # planewavecutoff in eV
calc.SetDensityCutoff(400)
calc.SetNumberOfBands(8)             # set the number of electronic bands
calc.SetSymmetryOn()                 # use symmetry to reduce the k-point set
calc.SetNetCDFFile('make_density_converge.nc')
calc.SetTxtFile('make_density_converge.text')

bulk.SetCalculator(calc)

calc.StayAliveOff()

energy = calc.GetPotentialEnergy()
fermilevel = calc.GetFermiLevel()

# setup the Harris calculation

# setup the line from the Gamme point (0,0,0) to the X point (0.5,0.5,0) 
N = 80
Nf = float(N) 
kpts = [(x/Nf,x/Nf,0) for x in range(N/2)]
harris = Dacapo()
harris.SetBZKPoints(kpts)          # set the k-points along one direction
harris.SetPlaneWaveCutoff(340)     # planewavecutoff in eV
harris.SetNumberOfBands(6)         # set the number of electronic bands
harris.SetSymmetryOn()             # use symmetry to reduce the k-point set
harris.SetNetCDFFile('calc_eigenvalues.nc')
harris.SetTxtFile('calc_eigenvalues.text')

bulk.SetCalculator(harris) 

harris.StayAliveOff()
harris.SetChargeMixing(False) 

# set the density found using the uniform k-point sampling
harris.SetDensityArray(calc.GetDensityArray())

energy1 = harris.GetPotentialEnergy()
kpoints     = harris.GetIBZKPoints()
eigenvalues = [harris.GetEigenvalues(kpt=kpt) for kpt in range(len(kpoints))]

# using the fermilevel from uniform k-point sampling
print fermilevel
print kpoints
print eigenvalues

</code>
===== 一回収斂させた電荷密度を使って計算する =====

<code python>
""" Make a bandstructure plot using the Harris functional. 
"""
import sys
from Dacapo import Dacapo
from ASE import Atom,ListOfAtoms
from ASE.Utilities.List import List
import Numeric as num
import os

bulk = Dacapo.ReadAtoms('density_converged.nc')
calc = bulk.GetCalculator()

energy = calc.GetPotentialEnergy()
fermilevel = calc.GetFermiLevel()

# setup the Harris calculation

# setup the line from the Gamme point (0,0,0) to the X point (0.5,0.5,0) 
N = 80
Nf = float(N) 
kpts = [(x/Nf,x/Nf,0) for x in range(N/2)]
harris = Dacapo()
harris.SetBZKPoints(kpts)          # set the k-points along one direction
harris.SetPlaneWaveCutoff(340)     # planewavecutoff in eV
harris.SetNumberOfBands(6)         # set the number of electronic bands
harris.SetSymmetryOn()             # use symmetry to reduce the k-point set
harris.SetNetCDFFile('calc_eigenvalues_only.nc')
harris.SetTxtFile('calc_eigenvalues_only.text')

bulk.SetCalculator(harris) 

harris.StayAliveOff()
harris.SetChargeMixing(False) 

# set the density found using the uniform k-point sampling
harris.SetDensityArray(calc.GetDensityArray())

energy1 = harris.GetPotentialEnergy()
kpoints     = harris.GetIBZKPoints()
eigenvalues = [harris.GetEigenvalues(kpt=kpt) for kpt in range(len(kpoints))]

# using the fermilevel from uniform k-point sampling
print fermilevel
print kpoints
print eigenvalues
</code>

<code python>

# -*- coding: iso-8859-15 -*-
# 日本語でコメントを書くために上の行が必要
from math import sqrt  # sqrt(x)で√xが計算できる
##################################################
ncfile = "a2.nc"  # このファイルは上書きされてしまうので
                        # 必ずコピーを取っておくこと

M = 10
Mf = float(M) 
wk1 = [(2*n/(3*Mf),n/(3*Mf),0) for n in range(M)]

from Dacapo import Dacapo
s2 = Dacapo.ReadAtoms(ncfile)
c2 = s2.GetCalculator()

c2.SetBZKPoints(wk1)    ### Band計算に必要
# c2.SetNetCDFFile("b1.nc")
c2.SetTxtFile("b2.text")
c2.StayAliveOff()         ### Band計算に必要
c2.SetChargeMixing(False) ### Band計算に必要 
c2.Calculate() ### Band計算に必要
</code>