Source code for aiida_crystal17.parsers.raw.parse_bases

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Copyright 2019 Chris Sewell
#
# This file is part of aiida-crystal17.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms and conditions
# of version 3 of the GNU Lesser General Public License.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
from collections import namedtuple
import copy

from aiida_crystal17.common.atoms import (
    ELECTRON_CONFIGURATIONS,
    GAUSSIAN_ORBITALS,
    SYMBOLS,
    SYMBOLS_R,
)


[docs]def parse_bsets_stdin(content, allow_comments=False, isolated=False): """parse basis sets from a crystal intput file Parameters ---------- content : str file content to parse allow_comments : bool if True, comments will be stripped before parsing isolated : bool if the basis sets are not within the crystal input file Returns ------- dict {'bs': {<atom_type>: [{'type': <type>, 'functions': [...]}, ...]}, 'ecp': {<atom_type>: [[...], ...]}} Raises ------ IOError if an error occurs during the parsing NotImplementedError if more than 2 basis sets / pseudopotentials are set for one atom type Notes ----- Standard Basis Set:: a_number_id n_shells # for each shell type shell_type n_functions charge scale_factor # if type=0, for n_functions exponent contraction_coefficient The atomic number Z is given by the remainder of the division of the conventional atomic number by 100 (2 max per species in positions not symmetry-related): - a_number_id < 200 > 1000: all electron basis set - a_number_id > 200 < 1000: valence electron basis set Valence-electron only calculations can be performed with the aid of effective core pseudo-potentials (ECP). The ECP input must be inserted into the basis set input of the atoms with conventional atomic number>200. Effective core pseudo-potentials (ECP) section (p. 75):: INPUT / HAYWLC / HAYWSC / BARTHE / DURAND # if INPUT insert effective_core_charge M M0 M1 M2 M3 M4 # insert M+M0+M1+M2+M3+M4 records a_kl C_kl n_kl """ gbasis = {} if not content: raise IOError("content is empty") comment_signals = "#/*<!" bs_sequence = {0: "S", 1: "SP", 2: "P", 3: "D", 4: "F", 5: "G"} bs_type = { 1: "STO-nG(nd) type ", # Pople standard STO-nG (Z=1-54); 2: "3(6)-21G(nd) type ", # Pople standard 3(6)-21G (Z=1-54(18)) + standard polarization functions } bs_notation = { 1: "n-21G outer valence shell", 2: "n-21G inner valence shell", 3: "3-21G core shell", 6: "6-21G core shell", } ps_keywords = { "INPUT": None, "HAYWLC": "Hay-Wadt large core", "HAYWSC": "Hay-Wadt small core", "BARTHE": "Durand-Barthelat", "DURAND": "Durand-Barthelat", } ps_sequence = ["W0", "P0", "P1", "P2", "P3", "P4"] in_basis_section = isolated in_pseudo, in_basis = False, False atomic_symbol = atomic_number_id = distrib = ps_indices_map = None for lineno, line in enumerate(content.splitlines(), 1): # ignore input until reaching the end of the geometry/optimisation section if line.startswith("END"): in_basis_section = True continue if not in_basis_section: continue # strip end of line comments if allow_comments: for s in comment_signals: pos = line.find(s) if pos != -1: line = line[:pos] break parts = line.split() if len(parts) == 0: continue if len(parts) == 1 and parts[0].upper() in ps_keywords: # start of pseudo if atomic_symbol is None: raise IOError( "line {}; reached pseudo before setting atom type".format(lineno) ) if not 200 < atomic_number_id < 999: raise IOError( "line {}; the basis contains an ecp, but the atomic_id {} is not in the range [200, 999]".format( lineno, atomic_number_id ) ) if parts[0] != "INPUT": gbasis[atomic_symbol]["ecp"].append(ps_keywords[parts[0].upper()]) continue # sanity check try: [float(p.replace("D", "E")) for p in parts] except ValueError: in_basis_section = False continue if len(parts) in [2, 3]: if parts[0] == "99" and parts[1] == "0": # end of basis set section break elif "." in parts[0] or "." in parts[1]: # exponent section for bs or ecp parts = [float(x.replace("D", "E")) for x in parts] if in_pseudo: # distribute exponents into ecp-types according to counter, that we now calculate if distrib in list(ps_indices_map.keys()): gbasis[atomic_symbol]["ecp"].append([ps_indices_map[distrib]]) gbasis[atomic_symbol]["ecp"][-1].append(tuple(parts)) distrib += 1 elif in_basis: # distribute exponents into orbitals according to counter, that we already defined gbasis[atomic_symbol]["bs"][-1]["functions"].append(tuple(parts)) else: # atom type definition section atomic_number_id = int(parts[0]) atomic_number = atomic_number_id % 100 if atomic_number == 0: atomic_symbol = "X" else: if atomic_number not in SYMBOLS: raise IOError( "line {}; atomic number not recognised: {}".format( lineno, atomic_number ) ) atomic_symbol = SYMBOLS[atomic_number] if atomic_symbol in gbasis: if atomic_symbol + "1" in gbasis: raise NotImplementedError( "line {}; More than two different basis sets for element {}".format( lineno, atomic_symbol ) ) atomic_symbol += "1" if 200 < atomic_number_id < 999: gbasis[atomic_symbol] = { "type": "valence-electron", "bs": [], "ecp": [], } else: gbasis[atomic_symbol] = {"type": "all-electron", "bs": []} elif len(parts) == 5: # orbital section gbasis[atomic_symbol]["bs"].append( { "type": bs_sequence[int(parts[1])], "functions": [], "charge": float(parts[3]), } ) parts = list(map(int, parts[0:3])) if parts[0] == 0: # insert from data given in input in_pseudo, in_basis = False, True elif parts[0] in bs_type: # pre-defined insert if parts[2] in bs_notation: gbasis[atomic_symbol]["bs"][-1]["functions"].append( bs_type[parts[0]] + bs_notation[parts[2]] ) else: gbasis[atomic_symbol]["bs"][-1]["functions"].append( bs_type[parts[0]] + "n=" + str(parts[2]) ) elif 6 <= len(parts) <= 7: # pseudo - INPUT section parts.pop(0) ps_indices = list(map(int, parts)) ps_indices_map = {} accum = 1 for c, n in enumerate(ps_indices): if n == 0: continue ps_indices_map[accum] = ps_sequence[c] accum += n distrib = 1 in_pseudo, in_basis = True, False for k, v in gbasis.items(): # sometimes no BS for host atom is printed when it is replaced by Xx: account for it if not len(v["bs"]) and k != "X" and "X" in gbasis: gbasis[k] = copy.deepcopy(gbasis["X"]) # NOTE: no GHOST deletion should be performed, since it breaks orbitals order for band structure plotting return gbasis
OrbitalResult = namedtuple( "OrbitalResult", ["electrons", "core_electrons", "number_ao", "orbital_types", "ao_indices"], )
[docs]def compute_orbitals(atoms, basis_sets): # type: (list, dict) -> OrbitalResult """compute data for all atomic orbitals in a structure, given elemental representations by crystal basis sets Parameters ---------- atoms : list[str] or list[int] list of atomic numbers or symbols which the structure comprises of basis_sets : dict[str, dict] basis set data, in the format returned from ``parse_bsets_stdin`` Returns ------- OrbitalResult """ total_electrons = total_core_electrons = total_aos = 0 aos_indices = {} orbital_types = [] for atom_index, atom in enumerate(atoms): try: electrons = int(atom) symbol = SYMBOLS[int(atom)] except (TypeError, ValueError): symbol = atom electrons = SYMBOLS_R[atom] if basis_sets[symbol]["type"] == "valence-electron": raise NotImplementedError("computing for bases with core pseudopotentials") outer_electrons = sum( [i for n, i in ELECTRON_CONFIGURATIONS[electrons]["outer"]] ) total_electrons += electrons total_core_electrons += electrons - outer_electrons type_count = {} for orbital in basis_sets[symbol]["bs"]: type_count.setdefault(orbital["type"], 0) type_count[orbital["type"]] += 1 for i in range(GAUSSIAN_ORBITALS[orbital["type"]]): total_aos += 1 if ( symbol, orbital["type"], type_count[orbital["type"]], ) not in orbital_types: orbital_types.append( (symbol, orbital["type"], type_count[orbital["type"]]) ) aos_indices[total_aos] = { "atom": atom_index, "element": symbol, "type": orbital["type"], "index": type_count[orbital["type"]], } return OrbitalResult( total_electrons, total_core_electrons, total_aos, orbital_types, aos_indices )