Properties: Density of States (DOSS)¶
The CryDossCalculation
can be used to run the properties
executable for DOSS calculations, from an existing fort.9
.
See also
Properties Workflow to run multiple properties calculations (and optional initial SCF).
!verdi plugin list aiida.calculations crystal17.doss
Inputs
code: required Code The Code to use for this job.
parameters: required Dict the input parameters to create the properties input file.
wf_folder: required FolderData, RemoteData, SinglefileData the folder containing the wavefunction fort.9 file
metadata: optional
Outputs
remote_folder: required RemoteData Input files necessary to run the process will be stored in this folder node ...
results: required Dict Summary Data extracted from the output file(s)
retrieved: required FolderData Files that are retrieved by the daemon will be stored in this node. By defa ...
arrays: optional ArrayData energies and DoS arrays
Exit codes
1: The process has failed with an unspecified error.
2: The process failed with legacy failure mode.
10: The process returned an invalid output.
11: The process did not register a required output.
200: The retrieved folder data node could not be accessed.
210: The main (stdout) output file was not found
211: The temporary retrieved folder was not found
300: An error was flagged trying to parse the crystal exec stdout file
350: The input file could not be read by crystal
351: Crystal could not find the required wavefunction file
352: Parser could not find the output isovalue (fort.25) file
353: Error parsing output isovalue (fort.25) file
400: The calculation stopped prematurely because it ran out of walltime.
401: The calculation stopped prematurely because it ran out of memory.
402: The calculation stopped prematurely because it ran out of virtual memory.
413: An error encountered usually during geometry optimisation
414: An error was encountered during an scf computation
415: An unknown error was encountered, causing the mpi to abort
499: The main crystal output file flagged an unhandled error
from aiida import engine, load_profile, orm, plugins
from aiida.tools.visualization import Graph
from aiida_crystal17.common import recursive_round
from aiida_crystal17.tests.utils import (
get_or_create_local_computer, get_or_create_code,
get_default_metadata)
from aiida_crystal17.tests import open_resource_binary
import matplotlib.pyplot as plt
%matplotlib inline
profile = load_profile()
profile.name
'test_crystal17'
Running a calculation¶
The DOSS Input Schema gives the allowed format of the input dictionary, for example:
computer = get_or_create_local_computer('work_directory', 'localhost')
code = get_or_create_code('crystal17.doss', computer, 'mock_properties17')
builder = code.get_builder()
builder.metadata = get_default_metadata()
builder.parameters = orm.Dict(dict={
"k_points": [18, 36],
"npoints": 100,
"band_minimum": -10,
"band_maximum": 10,
"band_units": "eV"
})
with open_resource_binary('doss', 'mgo_sto3g_scf', 'fort.9') as handle:
builder.wf_folder = orm.SinglefileData(handle)
result, calcnode = engine.run_get_node(builder)
!verdi process show {calcnode.pk}
Property Value
------------- ------------------------------------
type CalcJobNode
pk 884
uuid ff55e9d9-9a03-4375-b5d9-2353794f1c92
label
description
ctime 2019-09-24 14:04:37.274585+00:00
mtime 2019-09-24 14:04:51.167679+00:00
process state Finished
exit status 0
computer [1] localhost
Inputs PK Type
---------- ---- --------------
code 881 Code
parameters 882 Dict
wf_folder 883 SinglefileData
Outputs PK Type
------------- ---- ----------
arrays 888 ArrayData
remote_folder 885 RemoteData
results 887 Dict
retrieved 886 FolderData
graph = Graph(graph_attr={'size': "6,8!", "rankdir": "LR"})
graph.add_node(calcnode)
graph.add_incoming(calcnode, annotate_links="both")
graph.add_outgoing(calcnode, annotate_links="both")
graph.graphviz
Analysing the outputs¶
The outputs are:
results a dict of computation input and output parameters, parsed from the stdout file.
arrays a set of energy and DoS arrays for each point
recursive_round(calcnode.outputs.results.get_dict(), 1)
{'newk': {'k_points': [18, 18, 18],
'gilat_net': 36,
'n_kpoints_ibz': 195,
'n_kpoints_gilat': 1240},
'npts': 102,
'spin': False,
'units': {'energy': 'eV', 'conversion': 'CODATA2014'},
'errors': [],
'header': {'crystal_version': 17, 'crystal_subversion': '1.0.1'},
'warnings': [],
'wf_input': {'n_ao': 14,
'n_atoms': 2,
'k_points': [8, 8, 8],
'n_shells': 5,
'n_symops': 48,
'gilat_net': 8,
'n_core_el': 12,
'n_electrons': 20,
'energy_fermi': -4.0,
'energy_total': -7380.2,
'n_kpoints_ibz': 29,
'energy_kinetic': 7269.0},
'energy_max': 10.4,
'energy_min': -10.2,
'system_type': 'closed shell, insulating system',
'fermi_energy': -4.0,
'parser_class': 'CryDossParser',
'parser_errors': [],
'parser_version': '0.11.0',
'norbitals_total': 14,
'parser_exceptions': [],
'execution_time_seconds': 0}
plt.plot(
calcnode.outputs.arrays.get_array("energies"),
calcnode.outputs.arrays.get_array("total"))
plt.gca().set_xbound(-10, 0)
Computing Projections¶
Projections can be added per atom or per orbital set.
Note
A maximum of 15 projections are allowed per calculation.
Note
orm.Dict(dict={
"shrink_is": 18,
"shrink_isp": 36,
"npoints": 100,
"band_minimum": -10,
"band_maximum": 10,
"band_units": "eV",
"atomic_projections": [0, 1],
"orbital_projections": [[1, 2, 3]]
})
<Dict: uuid: c2ec5f24-bd4b-46fe-ba21-23b121ad3f6d (unstored)>
In order to create orbital sets,
it is possible to compute the nature of each orbital,
using the atomic structure and basis sets used to create the fort.9
:
from aiida_crystal17.tests import get_test_structure_and_symm
from aiida_crystal17.symmetry import print_structure
structure, _ = get_test_structure_and_symm('NiO_afm')
print_structure(structure)
StructureData Summary
Lattice
abc : 2.944 2.944 4.164
angles : 90.0 90.0 90.0
volume : 36.1
pbc : True True True
A : 2.944 0.0 0.0
B : 0.0 2.944 0.0
C : 0.0 0.0 4.164
Kind Symbols Position
---- ------- --------
Ni1 Ni 0.0 0.0 0.0
Ni2 Ni 1.472 1.472 2.082
O O 0.0 0.0 2.082
O O 1.472 1.472 0.0
basis_cls = plugins.DataFactory('crystal17.basisset')
basis_sets = basis_cls.get_basissets_from_structure(structure, 'sto3g')
basis_data = {k: v.get_data() for k, v in basis_sets.items()}
basis_data
{'Ni': {'type': 'all-electron',
'bs': [{'type': 'S', 'functions': ['STO-nG(nd) type 3-21G core shell']},
{'type': 'SP', 'functions': ['STO-nG(nd) type 3-21G core shell']},
{'type': 'SP', 'functions': ['STO-nG(nd) type 3-21G core shell']},
{'type': 'SP', 'functions': ['STO-nG(nd) type 3-21G core shell']},
{'type': 'D', 'functions': ['STO-nG(nd) type 3-21G core shell']}]},
'O': {'type': 'all-electron',
'bs': [{'type': 'S', 'functions': ['STO-nG(nd) type 3-21G core shell']},
{'type': 'SP', 'functions': ['STO-nG(nd) type 3-21G core shell']}]}}
from aiida_crystal17.parsers.raw.parse_bases import compute_orbitals
cresult = compute_orbitals(structure.get_ase().numbers, basis_data)
print("number of electrons: ", cresult.electrons)
print("number of core electrons: ", cresult.core_electrons)
cresult.ao_indices
number of electrons: 72
number of core electrons: 40
{1: {'atom': 0, 'element': 'Ni', 'type': 'S', 'index': 1},
2: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 1},
3: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 1},
4: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 1},
5: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 1},
6: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 2},
7: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 2},
8: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 2},
9: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 2},
10: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 3},
11: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 3},
12: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 3},
13: {'atom': 0, 'element': 'Ni', 'type': 'SP', 'index': 3},
14: {'atom': 0, 'element': 'Ni', 'type': 'D', 'index': 1},
15: {'atom': 0, 'element': 'Ni', 'type': 'D', 'index': 1},
16: {'atom': 0, 'element': 'Ni', 'type': 'D', 'index': 1},
17: {'atom': 0, 'element': 'Ni', 'type': 'D', 'index': 1},
18: {'atom': 0, 'element': 'Ni', 'type': 'D', 'index': 1},
19: {'atom': 1, 'element': 'Ni', 'type': 'S', 'index': 1},
20: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 1},
21: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 1},
22: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 1},
23: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 1},
24: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 2},
25: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 2},
26: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 2},
27: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 2},
28: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 3},
29: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 3},
30: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 3},
31: {'atom': 1, 'element': 'Ni', 'type': 'SP', 'index': 3},
32: {'atom': 1, 'element': 'Ni', 'type': 'D', 'index': 1},
33: {'atom': 1, 'element': 'Ni', 'type': 'D', 'index': 1},
34: {'atom': 1, 'element': 'Ni', 'type': 'D', 'index': 1},
35: {'atom': 1, 'element': 'Ni', 'type': 'D', 'index': 1},
36: {'atom': 1, 'element': 'Ni', 'type': 'D', 'index': 1},
37: {'atom': 2, 'element': 'O', 'type': 'S', 'index': 1},
38: {'atom': 2, 'element': 'O', 'type': 'SP', 'index': 1},
39: {'atom': 2, 'element': 'O', 'type': 'SP', 'index': 1},
40: {'atom': 2, 'element': 'O', 'type': 'SP', 'index': 1},
41: {'atom': 2, 'element': 'O', 'type': 'SP', 'index': 1},
42: {'atom': 3, 'element': 'O', 'type': 'S', 'index': 1},
43: {'atom': 3, 'element': 'O', 'type': 'SP', 'index': 1},
44: {'atom': 3, 'element': 'O', 'type': 'SP', 'index': 1},
45: {'atom': 3, 'element': 'O', 'type': 'SP', 'index': 1},
46: {'atom': 3, 'element': 'O', 'type': 'SP', 'index': 1}}
To observe DoS at the fermi level, these results can also be used to choose a sensible range of bands:
filled_bands = int(cresult.electrons / 2)
first_band = int(cresult.core_electrons / 2) + 1
last_band = min([first_band + 2 * (filled_bands - first_band), cresult.number_ao])
first_band, last_band
(21, 46)
O_sp_orbitals = [k for k, v in cresult.ao_indices.items()
if v["element"] == "O" and v["type"] == "SP"]
Ni_d_orbitals = [k for k, v in cresult.ao_indices.items()
if v["element"] == "Ni" and v["type"] == "D"]
computer = get_or_create_local_computer('work_directory', 'localhost')
code = get_or_create_code('crystal17.doss', computer, 'mock_properties17')
builder = code.get_builder()
builder.metadata = get_default_metadata()
builder.parameters = orm.Dict(dict={
"k_points": [18, 36],
"npoints": 1000,
"band_minimum": first_band,
"band_maximum": last_band,
"band_units": "bands",
"orbital_projections": [O_sp_orbitals, Ni_d_orbitals]
})
with open_resource_binary('doss', 'nio_sto3g_afm', 'fort.9') as handle:
builder.wf_folder = orm.SinglefileData(handle)
result2, calcnode2 = engine.run_get_node(builder)
calcnode2.outputs.arrays.attributes
{'array|energies': [1002],
'array|total_beta': [1002],
'array|total_alpha': [1002],
'array|projections_beta': [2, 1002],
'array|projections_alpha': [2, 1002]}
round(calcnode2.res.fermi_energy, 2)
6.89
plt.plot(
calcnode2.outputs.arrays.get_array("energies"),
calcnode2.outputs.arrays.get_array("total_alpha"),
label="Total $\\alpha$",
color="black", linestyle="dashed"
)
plt.plot(
calcnode2.outputs.arrays.get_array("energies"),
calcnode2.outputs.arrays.get_array("total_beta"),
label="Total $\\beta$",
color="black", linestyle="dashed"
)
plt.plot(
calcnode2.outputs.arrays.get_array("energies"),
calcnode2.outputs.arrays.get_array("projections_alpha")[0],
label="O(sp) $\\alpha$"
)
plt.plot(
calcnode2.outputs.arrays.get_array("energies"),
calcnode2.outputs.arrays.get_array("projections_alpha")[1],
label="Mg(d) $\\alpha$"
)
plt.plot(
calcnode2.outputs.arrays.get_array("energies"),
calcnode2.outputs.arrays.get_array("projections_beta")[0],
label="O(sp) $\\beta$"
)
plt.plot(
calcnode2.outputs.arrays.get_array("energies"),
calcnode2.outputs.arrays.get_array("projections_beta")[1],
label="Mg(d) $\\beta$"
)
plt.legend();
plt.gcf().set_size_inches(10, 8)