File: homo.py

package info (click to toggle)
gpaw 21.1.0-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 14,492 kB
  • sloc: python: 121,997; ansic: 14,138; sh: 1,125; csh: 139; makefile: 43
file content (61 lines) | stat: -rw-r--r-- 1,673 bytes parent folder | download | duplicates (2)
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
49
50
51
52
53
54
55
56
57
58
59
60
61
from ase.visualize import view
from ase.build import fcc111, add_adsorbate
from gpaw import GPAW
from gpaw.mixer import MixerSum
import gpaw.dscf as dscf

filename = 'homo'

c_mol = GPAW(nbands=9,
             h=0.2,
             xc='RPBE',
             kpts=(8, 6, 1),
             spinpol=True,
             convergence={'energy': 100,
                          'density': 100,
                          'eigenstates': 1.0e-9,
                          'bands': 'occupied'},
             txt='CO_homo.txt')

calc = GPAW(nbands=80,
            h=0.2,
            xc='RPBE',
            kpts=(8, 6, 1),
            eigensolver='cg',
            spinpol=True,
            mixer=MixerSum(nmaxold=5, beta=0.1, weight=100),
            convergence={'energy': 100,
                         'density': 100,
                         'eigenstates': 1.0e-7,
                         'bands': -10},
            txt=filename + '.txt')

# Import Slab with relaxed CO
slab = fcc111('Pt', size=(1, 2, 3), orthogonal=True)
add_adsorbate(slab, 'C', 2.0, 'ontop')
add_adsorbate(slab, 'O', 3.15, 'ontop')
slab.center(axis=2, vacuum=4.0)

view(slab)

molecule = slab.copy()

del molecule[:-2]

# Molecule
molecule.calc = c_mol
molecule.get_potential_energy()

# Homo wavefunction
wf_u = [kpt.psit_nG[4] for kpt in c_mol.wfs.kpt_u]

# Homo projector overlaps
mol = range(len(slab))[-2:]
p_uai = [dict([(mol[a], P_ni[4]) for a, P_ni in kpt.P_ani.items()])
         for kpt in c_mol.wfs.kpt_u]

# Slab with adsorbed molecule
slab.calc = calc
orbital = dscf.AEOrbital(calc, wf_u, p_uai, Estart=-100.0, Eend=0.0)
dscf.dscf_calculation(calc, [[-1.0, orbital, 1]], slab)
slab.get_potential_energy()