# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at https://mozilla.org/MPL/2.0/.
import os.path as osp
from numbers import Integral
import numpy as np
from sisl import Grid, Lattice
from sisl._internal import set_module
from sisl.messages import deprecate_argument, info
from sisl.unit.siesta import unit_convert
from .._help import grid_reduce_indices
from ..sile import SileError, add_sile, sile_raise_write
from .sile import SileCDFSiesta
__all__ = ['gridncSileSiesta']
Bohr2Ang = unit_convert('Bohr', 'Ang')
Ry2eV = unit_convert('Ry', 'eV')
@set_module("sisl.io.siesta")
class gridncSileSiesta(SileCDFSiesta):
""" NetCDF real-space grid file
The grid sile will automatically convert the units from Siesta units (Bohr, Ry) to sisl units (Ang, eV) provided the correct extension is present.
"""
[docs] def read_lattice(self):
""" Returns a Lattice object from a Siesta.grid.nc file
"""
cell = np.array(self._value('cell'), np.float64)
# Yes, this is ugly, I really should implement my unit-conversion tool
cell *= Bohr2Ang
cell.shape = (3, 3)
return Lattice(cell)
[docs] @deprecate_argument("sc", "lattice", "use lattice= instead of sc=", from_version="0.15")
def write_lattice(self, lattice):
""" Write a supercell to the grid.nc file """
sile_raise_write(self)
# Create initial dimensions
self._crt_dim(self, 'xyz', 3)
self._crt_dim(self, 'abc', 3)
v = self._crt_var(self, 'cell', 'f8', ('abc', 'xyz'))
v.info = 'Unit cell'
v.unit = 'Bohr'
v[:, :] = lattice.cell[:, :] / Bohr2Ang
[docs] def read_grid(self, index=0, name='gridfunc', *args, **kwargs):
""" Reads a grid in the current Siesta.grid.nc file
Enables the reading and processing of the grids created by Siesta
Parameters
----------
index : int or array_like, optional
the spin-index for retrieving one of the components. If a vector
is passed it refers to the fraction per indexed component. I.e.
``[0.5, 0.5]`` will return sum of half the first two components.
Default to the first component.
name : str, optional
the name for the grid-function (do not supply for standard Siesta output)
geometry: Geometry, optional
add the Geometry to the Grid
spin : optional
same as `index` argument. `spin` argument has precedence.
"""
# Default to *index* variable
index = kwargs.get("spin", index)
# Determine the name of this file
f = osp.basename(self.file)
# File names are made up of
# ElectrostaticPotential.grid.nc
# So the first one should be ElectrostaticPotential
try:
# <>.grid.nc
base = f.split('.')[-3]
except Exception:
base = 'None'
# Unit-conversion
BohrC2AngC = Bohr2Ang ** 3
unit = {'Rho': 1. / BohrC2AngC,
'DeltaRho': 1. / BohrC2AngC,
'RhoXC': 1. / BohrC2AngC,
'RhoInit': 1. / BohrC2AngC,
'Chlocal': 1. / BohrC2AngC,
'TotalCharge': 1. / BohrC2AngC,
'BaderCharge': 1. / BohrC2AngC,
'ElectrostaticPotential': Ry2eV,
'TotalPotential': Ry2eV,
'Vna': Ry2eV,
}.get(base, None)
# Fall-back
if unit is None:
unit = 1.
show_info = True
else:
show_info = False
# Swap as we swap back in the end
lattice = self.read_lattice().swapaxes(0, 2)
# Create the grid
nx = len(self._dimension('n1'))
ny = len(self._dimension('n2'))
nz = len(self._dimension('n3'))
if name is None:
v = self._variable('gridfunc')
else:
v = self._variable(name)
# Create the grid, Siesta uses periodic, always
grid = Grid([nz, ny, nx], bc=Grid.PERIODIC, lattice=lattice, dtype=v.dtype,
geometry=kwargs.get("geometry", None))
if v.ndim == 3:
grid.grid[:, :, :] = v[:, :, :] * unit
elif isinstance(index, Integral):
grid.grid[:, :, :] = v[index, :, :, :] * unit
else:
grid_reduce_indices(v, np.array(index) * unit, axis=0, out=grid.grid)
if show_info:
info(f"{self.__class__.__name__}.read_grid cannot determine the units of the grid. "
"The units may not be in sisl units.")
# Read the grid, we want the z-axis to be the fastest
# looping direction, hence x,y,z == 0,1,2
return grid.swapaxes(0, 2)
[docs] def write_grid(self, grid, spin=0, nspin=None, **kwargs):
""" Write a grid to the grid.nc file """
sile_raise_write(self)
# Default to *index* variable
spin = kwargs.get('index', spin)
self.write_lattice(grid.lattice)
if nspin is not None:
self._crt_dim(self, 'spin', nspin)
self._crt_dim(self, 'n1', grid.shape[0])
self._crt_dim(self, 'n2', grid.shape[1])
self._crt_dim(self, 'n3', grid.shape[2])
if nspin is None:
v = self._crt_var(self, "gridfunc", grid.dtype, ('n3', 'n2', 'n1'))
else:
v = self._crt_var(self, "gridfunc", grid.dtype, ('spin', 'n3', 'n2', 'n1'))
v.info = "Grid function"
if nspin is None:
v[:, :, :] = np.swapaxes(grid.grid, 0, 2)
else:
v[spin, :, :, :] = np.swapaxes(grid.grid, 0, 2)
add_sile('grid.nc', gridncSileSiesta)