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
../_images/calc_doss_10_0.svg

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)
../_images/calc_doss_14_0.png

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)
../_images/calc_doss_28_0.png