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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
from __future__ import print_function
import numpy as np
from ase.calculators.octopus import Octopus
from ase.collections import g2
from ase.build import bulk, graphene_nanoribbon
from ase.calculators.interfacechecker import check_interface
def calculate(name, system, **kwargs):
print('Calculate', name, system)
label = 'ink-%s' % name
kwargs0 = dict(stdout="'stdout.txt'",
FromScratch=True,
RestartWrite=False,
command='mpirun -np 4 octopus')
kwargs.update(**kwargs0)
calc = Octopus(label=label, **kwargs)
system.calc = calc
E = system.get_potential_energy()
eig = calc.get_eigenvalues()
check_interface(calc)
restartcalc = Octopus(label)
check_interface(restartcalc)
# Check reconstruction of Atoms object
new_atoms = restartcalc.get_atoms()
print('new')
print(new_atoms.positions)
calc2 = Octopus(label='ink-restart-%s' % name, **kwargs)
new_atoms.calc = calc2
E2 = new_atoms.get_potential_energy()
#print('energy', E, E2)
eig2 = calc2.get_eigenvalues()
eig_err = np.abs(eig - eig2).max()
e_err = abs(E - E2)
print('Restart E err', e_err)
print('Restart eig err', eig_err)
assert e_err < 1e-10
assert eig_err < 1e-10
return calc
if 1:
calc = calculate('H2O',
g2['H2O'],
OutputFormat='xcrysden',
Output='density + wfs + potential',
SCFCalculateDipole=True)
dipole = calc.get_dipole_moment()
E = calc.get_potential_energy()
print('dipole', dipole)
print('energy', E)
dipole_err = np.abs(dipole - [0., 0., -0.37]).max()
assert dipole_err < 0.02, dipole_err
energy_err = abs(-463.5944954 - E)
assert energy_err < 0.001, energy_err
if 1:
atoms = g2['O2']
atoms.center(vacuum=2.0)
calc = calculate('O2',
atoms,
BoxShape='parallelepiped',
SpinComponents='spin_polarized',
ExtraStates=2)
magmom = calc.get_magnetic_moment()
magmoms = calc.get_magnetic_moments()
print('magmom', magmom)
print('magmoms', magmoms)
if 1:
calc = calculate('Si',
bulk('Si', orthorhombic=True),
KPointsGrid=[[4, 4, 4]],
KPointsUseSymmetries=True,
Spacing=0.35)
eF = calc.get_fermi_level()
print('eF', eF)
if 0: # This calculation does not run will in Octopus
# We will do the "toothless" spin-polarised Si instead.
calc = calculate('Fe',
bulk('Fe', orthorhombic=True),
KPointsGrid=[[4, 4, 4]],
KPointsUseSymmetries=True,
ExtraStates=4,
Spacing=0.15,
SmearingFunction='fermi_dirac',
Smearing=0.1,
PseudoPotentialSet='sg15',
ExperimentalFeatures=True,
SpinComponents='spin_polarized')
eF = calc.get_fermi_level()
assert abs(eF - 5.33) < 1e-1
# XXXX octopus does not get magnetic state?
if 1:
calc = calculate('Si',
bulk('Si', orthorhombic=True),
KPointsGrid=[[4, 4, 4]],
SpinComponents='spin_polarized',
KPointsUseSymmetries=True,
Spacing=0.35)
eF = calc.get_fermi_level()
print('eF', eF)
if 0:
# Experimental feature: mixed periodicity. Let us not do this for now...
graphene = graphene_nanoribbon(2, 2, sheet=True)
graphene.positions = graphene.positions[:, [0, 2, 1]]
graphene.pbc = [1, 1, 0] # from 1, 0, 1
calc = calculate('graphene',
graphene,
KPointsGrid=[[2, 1, 2]],
KPointsUseSymmetries=True,
ExtraStates=4,
SmearingFunction='fermi_dirac',
Smearing=0.1)
|
