1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
from ase.build import molecule
from gpaw import GPAW
from gpaw import dscf
# Ground state calculation
calc = GPAW(nbands=8,
h=0.2,
xc='PBE',
spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -1})
CO = molecule('CO')
CO.center(vacuum=3)
CO.calc = calc
E_gs = CO.get_potential_energy()
# Obtain the pseudowavefunctions and projector overlaps of the
# state which is to be occupied. n=5,6 is the 2pix and 2piy orbitals
n = 5
molecule = [0, 1]
wf_u = [kpt.psit_nG[n] for kpt in calc.wfs.kpt_u]
p_uai = [dict([(molecule[a], P_ni[n]) for a, P_ni in kpt.P_ani.items()])
for kpt in calc.wfs.kpt_u]
# Excited state calculation
calc_es = GPAW(nbands=8,
h=0.2,
xc='PBE',
spinpol=True,
convergence={'energy': 100,
'density': 100,
'eigenstates': 1.0e-9,
'bands': -1})
CO.calc = calc_es
lumo = dscf.AEOrbital(calc_es, wf_u, p_uai)
# lumo = dscf.MolecularOrbital(calc, weights={0: [0, 0, 0, 1],
# 1: [0, 0, 0, -1]})
dscf.dscf_calculation(calc_es, [[1.0, lumo, 1]], CO)
E_es = CO.get_potential_energy()
print('Excitation energy: ', E_es - E_gs)
|
