====== Jacapoサンプル ======
このサンプルにはまだ実行できないものも含まれている。

===== Sample 1 =====

<code python>
#!/usr/bin/env python
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
co = Atoms([Atom('C',[0,0,0]),Atom('O',[1.2,0,0])],cell=(6.,6.,6.),pbc=True)
calc = Jacapo('sample1.nc',nbands=6,pw=350,ft=0.01,atoms=co)   
print co.get_potential_energy()
print co.get_forces()
</code>
===== Sample 2 =====
<code python>
#!/usr/bin/env python
from ase.calculators.jacapo import *
atoms = Jacapo.read_atoms('sample1.nc')
print atoms.get_potential_energy()
</code>

===== Sample 3 =====
<code python>
#!/usr/bin/env python
from ase import Atom, Atoms
from ase.calculators.jacapo import Jacapo
calc = Jacapo('sample1.nc')
ef = calc.get_ef() #get Fermi level
#it is common to normalize the energies so E_f is 0 eV
eigs = calc.get_eigenvalues()
occs  = calc.get_occupation_numbers()
print 'Eigenvalues = ',eigs 
print 'Occupancies = ',occs
#let's find the homo and lumo
occupied_state_energies = eigs[occs > 0.0]
unoccupied_state_energies = eigs[occs == 0.0]
homo = occupied_state_energies[-1]
lumo = unoccupied_state_energies[0]
print 'homo energy = %1.3f eV' % homo
print 'lumo energy = %1.3f eV' % lumo
#the fermi level is approximately half way between the homo and lumo
print 'E_f = %1.3f eV' % ef
print 'average homo/lumo = %1.3f eV' % ((homo+lumo)/2.)
print 'eigs', 'occs'
for i in range(len(eigs)):
  print eigs[i], occs[i]
</code>

===== Sample 4 =====
<code python>
#!/usr/bin/env python
from ase.dft.dos import *
from ase.calculators.jacapo import *
calc = Jacapo('sample1.nc')
dos = DOS(calc)
x = dos.get_energies()
y = dos.get_dos()
print 'Energy', 'DOS'
for i in range(len(x)):
  print x[i], y[i]

</code>

===== Sample 5 =====
<code python>
#!/usr/bin/env python
from ase.calculators.jacapo import *
calc = Jacapo('sample1.nc')
x,y,z,cd = calc.get_charge_density()

</code>

===== Sample 6 =====
<code python>
#!/usr/bin/env python
from ase.calculators.jacapo import *
calc = Jacapo('sample1.nc')
x,y,z,wf = calc.get_wave_function(band=5)

</code>

===== Sample 7 =====
<code python>
#!/usr/bin/env python
from ase import *
from ase.data import molecules
from ase.calculators.jacapo import *
### Setup calculators
benzene = Atoms('C6H6',positions = molecules.data['C6H6']['positions'],cell=[10,10,10],pbc=True)
benzene.center()
b1 = Jacapo('sample7_benzene.nc',nbands=18,pw=350,ft=0.01,atoms=benzene)  
chlorobenzene = Atoms('C6H6',positions = molecules.data['C6H6']['positions'],cell=[10,10,10],pbc=True)
chlorobenzene[11].set_symbol('Cl')
c1 = Jacapo('1.2.1.1-chlorobenzene.nc',nbands=18,pw=350,ft=0.01,atoms=chlorobenzene)   
x1,y1,z1,cd1 = b1.get_charge_density()
x2,y2,z2,cd2 = c1.get_charge_density()
cdiff = cd2 - cd1
print cdiff.min(),cdiff.max()

</code>

===== Sample 8 =====
<code python>
#!/usr/bin/env python
from ase import *
from ase.calculators.jacapo import *
np.set_printoptions(precision=3,suppress=True)

for d in [1.05, 1.1, 1.15, 1.2, 1.25]: #possible bond lengths
    co = Atoms([Atom('C',[0,0,0]),Atom('O',[d,0,0])],cell=(6,6,6),pbc=True)
    calc = Jacapo('sample8_%1.2f.nc' % d, pw=350, ft=0.01, atoms=co)
    print 'd = %1.2f ang' % d
    print 'energy = %f eV' % co.get_potential_energy()
    print 'forces = (eV/ang)\n', co.get_forces()
    print '' #blank line
    calc.strip() #removes large variables we don't need now.

</code>

===== Sample 9 =====
<code python>
#!/usr/bin/env python
from ase import *
from ase.calculators.jacapo import *
bondlengths = [1.05, 1.1, 1.15, 1.2, 1.25]
fnames = ['sample8_%1.2f.nc' % d for d in bondlengths]
atoms = [Jacapo.read_atoms(f) for f in fnames]
energies = [a.get_potential_energy() for a in atoms]
for i in range(len(bondlengths)):
  print bondlengths[i], energies[i]

</code>

===== Sample 10 =====
<code python>
#!/usr/bin/env python
from ase import *
from ase.calculators.jacapo import *
import numpy as np
co = Atoms([Atom('C',[0,0,0]),Atom('O',[1.2,0,0])],cell=(6,6,6),pbc=True)
calc = Jacapo('sample10.nc',pw=350,ft=0.01,atoms=co)  
qn = QuasiNewton(co, trajectory='sample10.traj')
qn.run()
pos = co.get_positions()
d = ((pos[0] - pos[1])**2).sum()**0.5
print 'Bondlength = %1.2f angstroms' % d
print co.get_forces()

</code>