ab_initio:jacapo_tutorial:lesson0
目次
Jacapoサンプル
このサンプルにはまだ実行できないものも含まれている。
Sample 1
#!/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()
Sample 2
#!/usr/bin/env python from ase.calculators.jacapo import * atoms = Jacapo.read_atoms('sample1.nc') print atoms.get_potential_energy()
Sample 3
#!/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]
Sample 4
#!/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]
Sample 5
#!/usr/bin/env python from ase.calculators.jacapo import * calc = Jacapo('sample1.nc') x,y,z,cd = calc.get_charge_density()
Sample 6
#!/usr/bin/env python from ase.calculators.jacapo import * calc = Jacapo('sample1.nc') x,y,z,wf = calc.get_wave_function(band=5)
Sample 7
#!/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()
Sample 8
#!/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.
Sample 9
#!/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]
Sample 10
#!/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()
ab_initio/jacapo_tutorial/lesson0.txt · 最終更新: 2022/08/23 13:34 by 127.0.0.1