from __future__ import print_function, division
import sys
from ctypes import POINTER, c_int64, c_float, c_char_p, create_string_buffer
import numpy as np
from scipy.sparse import csr_matrix
from pynao.m_libnao import libnao
from pynao.m_siesta_units import siesta_conv_coefficients
# interfacing with fortran subroutines
libnao.siesta_hsx_size.argtypes = (c_char_p, POINTER(c_int64), POINTER(c_int64),
POINTER(c_int64), POINTER(c_int64))
libnao.siesta_hsx_read.argtypes = (c_char_p, POINTER(c_int64), POINTER(c_float),
POINTER(c_int64), POINTER(c_int64),
POINTER(c_int64), POINTER(c_int64),
POINTER(c_int64))
# END of interfacing with fortran subroutines
[docs]def siesta_hsx_read(fname, force_gamma=None):
"""
Read the siesta.HSX file
"""
fname = create_string_buffer(fname.encode())
if force_gamma is None:
ft = c_int64(-1)
elif force_gamma:
ft = c_int64(1)
elif not force_gamma:
ft = c_int64(2)
bufsize, row_ptr_size, col_ind_size = c_int64(), c_int64(), c_int64()
libnao.siesta_hsx_size(fname, ft, bufsize, row_ptr_size, col_ind_size)
if bufsize.value<=0 or row_ptr_size.value <= 0 or col_ind_size.value <= 0:
return None
dat = np.empty(bufsize.value, dtype=np.float32)
dimensions = np.empty(4, dtype=np.int64)
row_ptr = np.empty(row_ptr_size.value, dtype=np.int64)
col_ind = np.empty(col_ind_size.value, dtype=np.int64)
libnao.siesta_hsx_read(fname, ft, dat.ctypes.data_as(POINTER(c_float)),
row_ptr.ctypes.data_as(POINTER(c_int64)), row_ptr_size,
col_ind.ctypes.data_as(POINTER(c_int64)), col_ind_size,
dimensions.ctypes.data_as(POINTER(c_int64)))
return dat, row_ptr, col_ind, dimensions
[docs]class siesta_hsx_c():
"""
Interface to the HSX data from Siesta calculations
"""
def __init__(self, **kw):
self.fname = fname = kw['fname'] if 'fname' in kw else 'siesta.HSX'
force_gamma = kw['force_gamma'] if 'force_gamma' in kw else None
dat, row_ptr, col_ind, dimensions = siesta_hsx_read(fname, force_gamma)
if dat is None or row_ptr is None or col_ind is None:
raise RuntimeError('file HSX not found '+ fname)
self.norbs, self.norbs_sc, self.nspin, self.nnz = dimensions
i = 0
self.is_gamma = (dat[i]>0); i=i+1;
self.nelec = int(dat[i]); i=i+1;
# The temperature from Siesta is given in Rydberg
self.telec = dat[i]*siesta_conv_coefficients["ry2ha"]; i=i+1;
self.h4 = np.reshape(dat[i:i+self.nnz*self.nspin], (self.nspin,self.nnz)); i=i+self.nnz*self.nspin;
self.s4 = dat[i:i+self.nnz]; i = i + self.nnz;
self.x4 = np.reshape(dat[i:i+self.nnz*3], (self.nnz,3)); i = i + self.nnz*3;
self.spin2h4_csr = []
for s in range(self.nspin):
self.spin2h4_csr.append(csr_matrix((self.h4[s,:], col_ind, row_ptr), dtype=np.float32))
self.s4_csr = csr_matrix((self.s4, col_ind, row_ptr), dtype=np.float32)
self.orb_sc2orb_uc=None
if(i < len(dat)):
if(self.is_gamma):
raise SystemError('i<len(dat) && gamma')
self.orb_sc2orb_uc = np.array(dat[i:i+self.norbs_sc]-1, dtype='int')
i = i + self.norbs_sc
if(i != len(dat)):
raise SystemError('i!=len(dat)')
[docs] def deallocate(self):
"""
Deallocate hsx variables
"""
del self.h4
del self.s4
del self.x4
del self.spin2h4_csr
del self.s4_csr