Skip to content

Commit

Permalink
ProjwfcParser: add support for non-collinear and spinorbit calculat…
Browse files Browse the repository at this point in the history
…ions (#546)

To support the new spin polarization modes, new subclasses of the data
tool `RealhydrogenOrbital`, which is based on the `OrbitalData` data
plugin, had to be created:

 * `NoncollinearHydrogenOrbital`: includes the spin projection `s_z`
   along the z-axis. Spin-up and spin-down states are treated as
   separate orbitals

 * `SpinorbitHydrogenOrbital`: includes the `total_angular_momentum`
   quantum number `j`, the `angular_momentum` quantum number `l` and the
   total `magnetic_number` `m_j`.
  • Loading branch information
lbotsch authored Sep 15, 2020
1 parent fe75c80 commit 031b92b
Show file tree
Hide file tree
Showing 31 changed files with 2,061 additions and 12 deletions.
90 changes: 78 additions & 12 deletions aiida_quantumespresso/parsers/projwfc.py
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,47 @@ def find_orbitals_from_statelines(out_info_dict):
:param out_info_dict: contains various technical internals useful in parsing
:return: orbitals, a list of orbitals suitable for setting ProjectionData
"""

# Format of statelines
# From PP/src/projwfc.f90: (since Oct. 8 2019)
#
# 1000 FORMAT (5x,"state #",i4,": atom ",i3," (",a3,"), wfc ",i2," (l=",i1)
# IF (lspinorb) THEN
# 1001 FORMAT (" j=",f3.1," m_j=",f4.1,")")
# ELSE IF (noncolin) THEN
# 1002 FORMAT (" m=",i2," s_z=",f4.1,")")
# ELSE
# 1003 FORMAT (" m=",i2,")")
# ENDIF
#
# Before:
# IF (lspinorb) THEN
# ...
# 1000 FORMAT (5x,"state #",i4,": atom ",i3," (",a3,"), wfc ",i2, &
# " (j=",f3.1," l=",i1," m_j=",f4.1,")")
# ELSE
# ...
# 1500 FORMAT (5x,"state #",i4,": atom ",i3," (",a3,"), wfc ",i2, &
# " (l=",i1," m=",i2," s_z=",f4.1,")")
# ENDIF

out_file = out_info_dict['out_file']
atomnum_re = re.compile(r'atom (.*?)\(')
element_re = re.compile(r'\((.*?)\)')
lnum_re = re.compile(r'l=(.*?)m=')
mnum_re = re.compile(r'm=(.*?)\)')
atomnum_re = re.compile(r'atom\s*([0-9]+?)[^0-9]')
element_re = re.compile(r'atom\s*[0-9]+\s*\(\s*([A-Za-z]+?)\s*\)')
if out_info_dict['spinorbit']:
# spinorbit
lnum_re = re.compile(r'l=\s*([0-9]+?)[^0-9]')
jnum_re = re.compile(r'j=\s*([0-9.]+?)[^0-9.]')
mjnum_re = re.compile(r'm_j=\s*([-0-9.]+?)[^-0-9.]')
elif not out_info_dict['collinear']:
# non-collinear
lnum_re = re.compile(r'l=\s*([0-9]+?)[^0-9]')
mnum_re = re.compile(r'm=\s*([-0-9]+?)[^-0-9]')
sznum_re = re.compile(r's_z=\s*([-0-9.]*?)[^-0-9.]')
else:
# collinear / no spin
lnum_re = re.compile(r'l=\s*([0-9]+?)[^0-9]')
mnum_re = re.compile(r'm=\s*([-0-9]+?)[^-0-9]')
wfc_lines = out_info_dict['wfc_lines']
state_lines = [out_file[wfc_line] for wfc_line in wfc_lines]
state_dicts = []
Expand All @@ -35,8 +71,14 @@ def find_orbitals_from_statelines(out_info_dict):
state_dict['atomnum'] -= 1 # to keep with orbital indexing
state_dict['kind_name'] = element_re.findall(state_line)[0].strip()
state_dict['angular_momentum'] = int(lnum_re.findall(state_line)[0])
state_dict['magnetic_number'] = int(mnum_re.findall(state_line)[0])
state_dict['magnetic_number'] -= 1 # to keep with orbital indexing
if out_info_dict['spinorbit']:
state_dict['total_angular_momentum'] = float(jnum_re.findall(state_line)[0])
state_dict['magnetic_number'] = float(mjnum_re.findall(state_line)[0])
else:
if not out_info_dict['collinear']:
state_dict['spin'] = float(sznum_re.findall(state_line)[0])
state_dict['magnetic_number'] = int(mnum_re.findall(state_line)[0])
state_dict['magnetic_number'] -= 1 # to keep with orbital indexing
except ValueError:
raise QEOutputParsingError('State lines are not formatted in a standard way.')
state_dicts.append(state_dict)
Expand All @@ -61,9 +103,14 @@ def find_orbitals_from_statelines(out_info_dict):

# here we set the resulting state_dicts to a new set of orbitals
orbitals = []
RealhydrogenOrbital = OrbitalFactory('realhydrogen')
if out_info_dict['spinorbit']:
OrbitalCls = OrbitalFactory('spinorbithydrogen')
elif not out_info_dict['collinear']:
OrbitalCls = OrbitalFactory('noncollinearhydrogen')
else:
OrbitalCls = OrbitalFactory('realhydrogen')
for state_dict in state_dicts:
orbitals.append(RealhydrogenOrbital(**state_dict))
orbitals.append(OrbitalCls(**state_dict))

return orbitals

Expand Down Expand Up @@ -190,6 +237,8 @@ def spin_dependent_pdos_subparser(out_info_dict):
:return: (pdos_name, pdos_array) tuples for all the specific pdos
"""
spin = out_info_dict['spin']
collinear = out_info_dict.get('collinear', True)
spinorbit = out_info_dict.get('spinorbit', False)
spin_down = out_info_dict['spin_down']
pdos_atm_array_dict = out_info_dict['pdos_atm_array_dict']
if spin:
Expand All @@ -210,8 +259,16 @@ def spin_dependent_pdos_subparser(out_info_dict):
# both are produced in the same order (thus the sorted file_names)
for name in pdos_file_names:
this_array = pdos_atm_array_dict[name]
for i in range(fa, np.shape(this_array)[1], mf):
out_arrays.append(this_array[:, i])
if not collinear and not spinorbit:
# In the non-collinear, non-spinorbit case, the "up"-spin orbitals
# come first, followed by all "down" orbitals
for i in range(3, np.shape(this_array)[1], 2):
out_arrays.append(this_array[:, i])
for i in range(4, np.shape(this_array)[1], 2):
out_arrays.append(this_array[:, i])
else:
for i in range(fa, np.shape(this_array)[1], mf):
out_arrays.append(this_array[:, i])

return out_arrays

Expand Down Expand Up @@ -262,7 +319,7 @@ def parse(self, **kwargs):
pdostot_filename = fnmatch.filter(out_filenames, '*pdos_tot*')[0]
with out_folder.open(pdostot_filename, 'r') as pdostot_file:
# Columns: Energy(eV), Ldos, Pdos
pdostot_array = np.genfromtxt(pdostot_file)
pdostot_array = np.atleast_2d(np.genfromtxt(pdostot_file))
energy = pdostot_array[:, 0]
dos = pdostot_array[:, 1]
except (OSError, KeyError):
Expand All @@ -273,7 +330,7 @@ def parse(self, **kwargs):
pdos_atm_array_dict = {}
for name in pdos_atm_filenames:
with out_folder.open(name, 'r') as pdosatm_file:
pdos_atm_array_dict[name] = np.genfromtxt(pdosatm_file)
pdos_atm_array_dict[name] = np.atleast_2d(np.genfromtxt(pdosatm_file))

# finding the bands and projections
# we create a dictionary the progressively accumulates more info
Expand Down Expand Up @@ -353,8 +410,17 @@ def _parse_bands_and_projections(self, out_info_dict):
spin = True
else:
spin = False
out_info_dict['spinorbit'] = parent_param.get_dict().get('spin_orbit_calculation', False)
out_info_dict['collinear'] = not parent_param.get_dict().get('non_colinear_calculation', False)
if not out_info_dict['collinear']:
# Sanity check
if nspin != 4:
raise QEOutputParsingError('The calculation is non-collinear, but nspin is not set to 4!')
spin = False
except KeyError:
spin = False
out_info_dict['spinorbit'] = False
out_info_dict['collinear'] = True
out_info_dict['spin'] = spin

# changes k-numbers to match spin
Expand Down
Empty file.
Empty file.
65 changes: 65 additions & 0 deletions aiida_quantumespresso/tools/data/orbital/noncollinearhydrogen.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@
# -*- coding: utf-8 -*-
"""A module defining hydrogen-like orbitals with non-collinear spin component."""

from aiida.common.exceptions import ValidationError

from aiida.tools.data.orbital.realhydrogen import RealhydrogenOrbital


def validate_spin(value):
"""Validate the value of the spin projection s_z."""
if not isinstance(value, (int, float)):
raise ValidationError('spin projection (s_z) must be float')

return value


class NoncollinearHydrogenOrbital(RealhydrogenOrbital):
"""Orbitals for hydrogen, with non-collinear spin.
The class largely follows the conventions used by wannier90 Object to
handle the generation of real hydrogen orbitals and their hybrids, has
methods for producing s, p, d, f, and sp, sp2, sp3, sp3d, sp3d2 hybrids.
This method does not deal with the cannonical hydrogen orbitals which
contain imaginary components. The orbitals described here are chiefly
concerned with the symmetric aspects of the oribitals without the context
of space. Therefore diffusitivity, position and atomic labels should be
handled in the OrbitalData class. Following the notation of table 3.1, 3.2
of Wannier90 user guide http://www.wannier.org/doc/user_guide.pdf A brief
description of what is meant by each of these labels: :param radial_nodes:
the number of radial nodes (or inflections) if no radial nodes are
supplied, defaults to 0 :param angular_momentum: Angular quantum number,
using real orbitals :param magnetic_number: Magnetic quantum number, using
real orbitals :param spin: spin z-projection s_z The conventions regarding
L and M correpsond to those used in wannier90 for all L greater than 0 the
orbital is not hyrbridized see table 3.1 and for L less than 0 the orbital
is hybridized see table 3.2. M then indexes all the possible orbitals from
0 to 2L for L > 0 and from 0 to (-L) for L < 0.
"""

_base_fields_required = RealhydrogenOrbital._base_fields_required

_base_fields_optional = tuple(
list(filter(lambda x: x[0] != 'spin', RealhydrogenOrbital._base_fields_optional)) + [
('spin', validate_spin, 0.0),
]
)

def __str__(self):
"""Printable representation of the orbital."""
orb_dict = self.get_orbital_dict()
try:
orb_name = self.get_name_from_quantum_numbers(
orb_dict['angular_momentum'], magnetic_number=orb_dict['magnetic_number']
)
position_string = '{:.4f},{:.4f},{:.4f}'.format(
orb_dict['position'][0], orb_dict['position'][1], orb_dict['position'][2]
)
out_string = 'r{} {} (s_z={}) orbital {} @ {}'.format(
orb_dict['radial_nodes'], orb_name, orb_dict['spin'],
"for kind '{}'".format(orb_dict['kind_name']) if orb_dict['kind_name'] else '', position_string
)
except KeyError:
# Should not happen, but we want it not to crash in __str__
out_string = '(not all parameters properly set)'
return out_string
143 changes: 143 additions & 0 deletions aiida_quantumespresso/tools/data/orbital/spinorbithydrogen.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,143 @@
# -*- coding: utf-8 -*-
"""A module defining hydrogen-like orbitals.
The orbitals are defined in the simultaneous eigen-basis of H, J^2, L^2, S^2
and J_z operators, useful for projecting densities with spin-orbit coupling.
"""

from aiida.common.exceptions import ValidationError

from aiida.tools.data.orbital.orbital import Orbital


def validate_j(value):
"""Validate the value of the total angular momentum."""
if not isinstance(value, (int, float)):
raise ValidationError('total angular momentum (j) must be float')

if not value % 0.5 == 0.0:
raise ValidationError('total angular momentum (j) must be a half integer')

if value < 0.5 or value > 3.5:
raise ValidationError('total angular momentum (j) must be between 0.5 and 3.5')

return value


def validate_l(value):
"""Validate the value of the angular momentum."""
if not isinstance(value, int):
raise ValidationError('angular momentum (l) must be integer')

if value < 0 or value > 3:
raise ValidationError('angular momentum (l) must be between 0 and 3')

return value


def validate_mj(value):
"""Validate the value of the magnetic number."""
if not isinstance(value, (int, float)):
raise ValidationError('magnetic number (m_j) must be float')

# without knowing j we cannot validate the value of m_j. We will call an additional function
# in the validate function

return value


def validate_kind_name(value):
"""Validate the value of the kind_name."""
if value is not None and not isinstance(value, str):
raise ValidationError('kind_name must be a string')

return value


def validate_n(value):
"""Validate the value of the number of radial nodes."""
if not isinstance(value, int):
raise ValidationError('number of radial nodes (n) must be integer')

if value < 0 or value > 2:
raise ValidationError('number of radial nodes (n) must be between 0 and 2')

return value


class SpinorbitHydrogenOrbital(Orbital):
"""Orbitals for hydrogen with spin-orbit interaction.
The orbital is defined in the common basis of H, J^2, L^2, S^2, J_z
operators, indexed by the quantum numbers n, j, l, j_z.
A brief description of what is meant by each of these labels:
:param radial_nodes: the number of radial nodes (or inflections) if no radial
nodes are supplied, defaults to 0
:param angular_momentum: Angular quantum number l
:param total_angular_momentum: Total angular momentum number j
:param magnetic_number: Magnetic quantum number m_j
The total angular momentum quantum number j takes values `|l-s| <= j <=
l+s` in integer steps and the magnetic number takes values from -j to +j in
integer steps.
"""

_base_fields_required = tuple(
list(Orbital._base_fields_required) +
[('total_angular_momentum',
validate_j), ('angular_momentum', validate_l), ('magnetic_number', validate_mj), ('radial_nodes', validate_n)]
)

_base_fields_optional = tuple(list(Orbital._base_fields_optional) + [
('kind_name', validate_kind_name, None),
])

def __str__(self):
"""Printable representation of the orbital."""
orb_dict = self.get_orbital_dict()
try:
orb_name = 'j={},l={},m_j={}'.format(
orb_dict['total_angular_momentum'], orb_dict['angular_momentum'], orb_dict['magnetic_number']
)
position_string = '{:.4f},{:.4f},{:.4f}'.format(
orb_dict['position'][0], orb_dict['position'][1], orb_dict['position'][2]
)
out_string = 'r{} {} orbital {} @ {}'.format(
orb_dict['radial_nodes'], orb_name,
"for kind '{}'".format(orb_dict['kind_name']) if orb_dict['kind_name'] else '', position_string
)
except KeyError:
# Should not happen, but we want it not to crash in __str__
out_string = '(not all parameters properly set)'
return out_string

def _validate_keys(self, input_dict):
"""Validate the keys otherwise raise ValidationError.
Does basic validation from the parent followed by validations for the
quantum numbers. Raises exceptions should the input_dict fail the
valiation or if it contains any unsupported keywords. :param
input_dict: the dictionary of keys to be validated :return
validated_dict: a validated dictionary
"""
validated_dict = super()._validate_keys(input_dict)

# Validate m knowing the value of l
total_angular_momentum = validated_dict['total_angular_momentum'] # j quantum number, must be there
angular_momentum = validated_dict['angular_momentum'] # l quantum number, must be there
accepted_range = [abs(angular_momentum - 0.5), angular_momentum + 0.5]
if total_angular_momentum < min(accepted_range) or total_angular_momentum > max(accepted_range):
raise ValidationError(
'the total angular momentum must be in the range [{}, {}]'.format(
min(accepted_range), max(accepted_range)
)
)
magnetic_number = validated_dict['magnetic_number'] # m quantum number, must be there
accepted_range = [-total_angular_momentum, total_angular_momentum]
if magnetic_number < min(accepted_range) or magnetic_number > max(accepted_range):
raise ValidationError(
'the magnetic number must be in the range [{}, {}]'.format(min(accepted_range), max(accepted_range))
)

return validated_dict
4 changes: 4 additions & 0 deletions setup.json
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,10 @@
"aiida.tools.calculations": [
"quantumespresso.pw = aiida_quantumespresso.tools.calculations.pw:PwCalculationTools"
],
"aiida.tools.data.orbitals": [
"spinorbithydrogen = aiida_quantumespresso.tools.data.orbital.spinorbithydrogen:SpinorbitHydrogenOrbital",
"noncollinearhydrogen = aiida_quantumespresso.tools.data.orbital.noncollinearhydrogen:NoncollinearHydrogenOrbital"
],
"aiida.tools.dbexporters.tcod_plugins": [
"quantumespresso.cp = aiida_quantumespresso.tools.dbexporters.tcod_plugins.cp:CpTcodtranslator",
"quantumespresso.pw = aiida_quantumespresso.tools.dbexporters.tcod_plugins.pw:PwTcodtranslator"
Expand Down
15 changes: 15 additions & 0 deletions tests/parsers/fixtures/projwfc/noncollinear/aiida.in
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
&PROJWFC
degauss = 2.0000000000d-02
deltae = 1.0000000000d-01
emax = 7.9029595890d+00
emin = 5.9029595890d+00
kresolveddos = .false.
lbinary_data = .false.
lsym = .true.
lwrite_overlaps = .false.
ngauss = 0
outdir = './out/'
plotboxes = .false.
prefix = 'aiida'
tdosinboxes = .false.
/
Loading

0 comments on commit 031b92b

Please sign in to comment.