""" 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

""" 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

# -*- 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計算に必要