import numpy as np
import scipy.sparse as sparse
import numba as nb
from pynao.m_vxc_sparse import vxc_csr_spin1, vxc_csr_spin2
[docs]def kernel_initialization(self, **kw):
"""
Load the kernel from file or compute it
Inputs:
self: scf (scf.py) or tddft_iter (tddft_iter.py) class
the contructor that needs kernel initialization. Should contains all
the necessary field to initialize the kernel
kernel_format: string
the storage format of the kernel, can be
* pack: store the kernel as a 1D array (lower part of the matrix)
* sparse: store the lower part of the kernel in a sparse format.
load_kernel: bool
Load the kernel from file (pack format only)
"""
if self.verbosity > 2:
print("{}\t====> kernel initialization".format(__name__))
if "xc_code" not in kw.keys():
raise ValueError("xc_code not in kw")
kernel_format = kw['kernel_format'] if 'kernel_format' in kw else "pack"
threshold = kw['kernel_threshold'] if 'kernel_threshold' in kw else 1.0e-14
load = kw["load_kernel"] if "load_kernel" in kw else False
xc_code = kw["xc_code"]
xc = xc_code.split(',')[0]
if load:
if xc != 'RPA' and self.nspin != 1:
raise RuntimeError('not sure it would work')
if kernel_format != "pack":
raise ValueError("only pack kernel format can be loaded")
params = {}
keys = ["kernel_fname", "kernel_file_format", "kernel_path_hdf5"]
for key in keys:
if key in kw.keys():
params[key] = kw[key]
kernel, kernel_dim = load_kernel_pack(self.nprod, self.dtype, **params)
else:
if kernel_format == "pack":
# lower triangular
kernel, kernel_dim = self.pb.comp_coulomb_pack(dtype=self.dtype)
elif kernel_format == "sparse":
kernel, kernel_dim = self.pb.comp_coulomb_sparse(dtype=self.dtype,
threshold=threshold)
else:
raise ValueError("unknow kernel_format {}".format(kernel_format))
assert self.nprod == kernel_dim, "{} {}".format(self.nprod, kernel_dim)
if self.verbosity > 2:
print("{}\t====> Hartree kernel calculated".format(__name__))
# List of POINTERS !!!
# of kernel [[(up,up), (up,dw)], [(dw,up), (dw,dw)]] TAKE CARE!!!
if xc == 'RPA' or xc == 'HF':
if self.nspin == 1:
ss2kernel = [[kernel]]
elif self.nspin == 2:
ss2kernel = [[kernel, kernel],
[kernel, kernel]]
elif xc == 'LDA' or xc == 'GGA':
kernel, ss2kernel = get_VXC_kernel(self, self.nspin, kernel, **kw)
else:
print(' xc_code', xc_code, xc, xc_code.split(','))
raise RuntimeError('unkn xc_code')
if self.verbosity > 2 and kernel_format == "sparse":
print("{}\t====> kernel size: {}; nnz: {}".format(__name__,
kernel.shape[0]*kernel.shape[1],
kernel.nnz))
if kernel_format == "sparse":
kernel_diag = kernel.diagonal()
ss2kernel_diag = []
if self.nspin == 2:
for spin1 in range(self.nspin):
ss2kernel_diag.append([])
for spin2 in range(self.nspin):
ss2kernel_diag[spin1].append(ss2kernel[spin1][spin2].diagonal())
else:
kernel_diag = None
ss2kernel_diag = None
return kernel_dim, kernel, kernel_diag, ss2kernel, ss2kernel_diag
[docs]def get_VXC_kernel(self, Nspin, kernel, **kw):
"""
Get the exchange-correlation part of the kernel and add it to Hartree kernel
Input parameters:
self: scf (scf.py) or tddft_iter (tddft_iter.py) class
the contructor that needs kernel initialization. Should contains all
the necessary field to initialize the kernel
Nspin: int
Spin number
kernel: 1D np.array or sparse matrix
The Hartree kernel
kernel_format: string
the storage format of the kernel, can be
* pack: store the kernel as a 1D array (lower part of the matrix)
* sparse: store the lower part of the kernel in a sparse format.
kernel_threshold: float
Threshold used for sparse kernel in order to drop small matrix element
and increase the sparsity of the matrix
"""
from pynao.m_vxc_pack import vxc_pack
from pynao.m_vxc_sparse import vxc_csr_spin1, vxc_csr_spin2
if self.verbosity > 2:
print("{}\t====> get_VXC_kernel".format(__name__))
kernel_format = kw['kernel_format'] if 'kernel_format' in kw else "pack"
threshold = kw['kernel_threshold'] if 'kernel_threshold' in kw else 1.0e-14
if Nspin == 1:
ss2kernel = [[None]]
if kernel_format == "pack":
vxc_pack(self, deriv=2, ao_log=self.pb.prod_log, kernel=kernel, **kw)
elif kernel_format == "sparse":
kernel += vxc_csr_spin1(self, deriv=2, ao_log=self.pb.prod_log,
threshold=threshold, **kw)
else:
raise ValueError("unknow kernel_format {}".format(kernel_format))
elif Nspin == 2:
if kernel_format == "pack":
kkk = vxc_pack(self, deriv=2, ao_log=self.pb.prod_log, **kw) + kernel
ss2kernel = [[kkk[0], kkk[1]],
[kkk[1], kkk[2]]]
elif kernel_format == "sparse":
kkk = vxc_csr_spin2(self, deriv=2, ao_log=self.pb.prod_log, threshold=threshold, **kw)
ss2kernel = [[kernel + kkk[0], kernel + kkk[1]],
[kernel + kkk[1], kernel + kkk[2]]]
else:
raise ValueError("unknow kernel_format {}".format(kernel_format))
return kernel, ss2kernel
[docs]@nb.jit(nopython=True)
def get_diagonal_triangle_matrix(arr, n):
diag = np.zeros((n))
idx = 0
for i in range(n):
for j in range(i+1):
if i == j:
diag[i] = arr[idx]
idx += 1
return diag
[docs]def load_kernel_pack(nprod, dtype, kernel_fname="", kernel_file_format="npy",
kernel_path_hdf5=""):
"""
Loads the kernel from file in a pack format and initializes field
Inputs parameters:
nprod: int
kernel shape is (nprod, nprod)
dtype: numpy data type
data type of the saved kernel, should be np.float32 or np.float64
kernel_fname: string
file name of the kernel
kernel_file_format: string
format of the file: txt, npy or hdf5
kernel_path_hdf5: string
where to find the kernel in the hdf5 file
Outputs parameters:
kernel: 1D np.array, float
the loaded kernel in pack format, size nprod*(nprod + 1)//2
kernel_dim: tuple, int
the dimension (shape) of the unpacked kernel, (nprod, nprod)
"""
if kernel_file_format == "npy":
kernel = dtype(np.load(kernel_fname))
elif kernel_file_format == "txt":
kernel = np.loadtxt(kernel_fname, dtype=dtype)
elif kernel_file_format == "hdf5":
import h5py
if not kernel_path_hdf5:
raise ValueError("kernel_path_hdf5 not set while trying to read kernel from hdf5 file.")
kernel = h5py.File(kernel_fname, "r")[kernel_path_hdf5][()]
else:
raise ValueError("Wrong format for loading kernel, must be: npy, txt or hdf5, got " \
+ kernel_file_format)
if len(kernel.shape) > 1:
raise ValueError("The kernel must be saved in packed format in order to be loaded!")
assert nprod*(nprod+1)//2 == kernel.size, \
"wrong size for loaded kernel: %r %r "%(nprod*(nprod+1)//2, kernel.size)
kernel_dim = nprod
return kernel, kernel_dim
[docs]def apply_kernel_nspin1_pack(spmv, nprod, kernel, dn, dtype=np.float64):
"""
Perform the matrix vector multiplication for the pack kernel with chi0_mv
in the iterative procedure
Input parameters:
spmv: blas function from scipy.linalg.blas
can be blas.sspmv for float or blas.dspmv for double
nprod: int
the dimension of the kernel matrix, i.e., (nprod, nprod)
kernel: 1D np.array, float or double
the kernel in pack format
dn: 1D np.array, complex
the resulting vector from chi0_mv
dtype: numpy data type
the data type of the kernel, np.float32 or np.float64
Output parameters:
vcre: 1D np.array, float
the real part of the matrix vector multiplication
vcim: 1D np.array, float
the imaginary part of the matrix vector multiplication
"""
daux = np.zeros(nprod, dtype=dtype)
vcreim = []
for dnreim in [dn.real, dn.imag]:
daux[:] = np.require(dnreim, dtype=dtype, requirements=["A","O"])
vcreim.append(spmv(nprod, 1.0, kernel, daux))
return vcreim[0], vcreim[1]
[docs]def apply_kernel_nspin2_pack(spmv, nprod, nspin, ss2kernel, dn, dtype=np.float64):
"""
Perform the matrix vector multiplication for the pack kernel with chi0_mv
in the iterative procedure for polarized spin calculation
Input parameters:
spmv: blas function from scipy.linalg.blas
can be blas.sspmv for float or blas.dspmv for double
nprod: int
the dimension of the kernel matrix, i.e., (nprod, nprod)
nspin: int
the number of spin
ss2kernel: list
the spin dependent kernel in pack format
ss2kernel = [[kkk[0], kkk[1]],
[kkk[1], kkk[2]]
dn: 1D np.array, complex
the resulting vector from chi0_mv
dtype: numpy data type
the data type of the kernel, np.float32 or np.float64
Output parameters:
vcre: 2D np.array, float or double
the real part of the matrix vector multiplication per spin
vcim: 2D np.array, float or double
the imaginary part of the matrix vector multiplication per spin
"""
vcre = np.zeros((2, nspin, nprod), dtype=dtype)
daux = np.zeros((nprod), dtype=dtype)
s2dn = dn.reshape((nspin, nprod))
for spin1 in range(nspin):
for spin2 in range(nspin):
for ireim, sreim in enumerate(('real', 'imag')):
daux[:] = np.require(getattr(s2dn[spin2], sreim), dtype=dtype,
requirements=["A","O"])
vcre[ireim, spin1] += spmv(nprod, 1.0,
ss2kernel[spin1][spin2], daux)
return vcre[0].reshape(-1), vcre[1].reshape(-1)
[docs]def apply_kernel_nspin1_sparse(nprod, kernel, kernel_diag, dn, dtype=np.float64):
"""
Perform the matrix vector multiplication for the sparse kernel with chi0_mv
in the iterative procedure
Input parameters:
nprod: int
the dimension of the kernel matrix, i.e., (nprod, nprod)
kernel: scipy csr_matrix, flat or double
the kernel in sparse format
dn: 1D np.array, complex
the resulting vector from chi0_mv
dtype: numpy data type
the data type of the kernel, np.float32 or np.float64
Output parameters:
vcre: 1D np.array, float or double
the real part of the matrix vector multiplication
vcim: 1D np.array, float or double
the imaginary part of the matrix vector multiplication
"""
daux = np.zeros(nprod, dtype=dtype)
vcreim = []
for dnreim in [dn.real, dn.imag]:
daux[:] = np.require(dnreim, dtype=dtype, requirements=["A","O"])
# Sparse triangular kermel matrix
# The matrix vector multiplication Adense.dot(x) writes as
#
# Asparse_tri.dot(x) + Asparse_tri.T.dot(x) - Asparse_tri.diagonal()*x
vcreim.append(kernel.dot(daux) + kernel.T.dot(daux) - kernel_diag*daux)
return vcreim[0], vcreim[1]
[docs]def apply_kernel_nspin2_sparse(nprod, nspin, ss2kernel, ss2kernel_diag, dn,
dtype=np.float64):
"""
Perform the matrix vector multiplication for the sparse kernel with chi0_mv
in the iterative procedure for polarized spin calculation
Input parameters:
nprod: int
the dimension of the kernel matrix, i.e., (nprod, nprod)
nspin: int
the number of spin
ss2kernel: list
the spin dependent kernel in sparse csr_matrix format
ss2kernel = [[kkk[0], kkk[1]],
[kkk[1], kkk[2]]
dn: 1D np.array, complex
the resulting vector from chi0_mv
dtype: numpy data type
the data type of the kernel, np.float32 or np.float64
Output parameters:
vcre: 2D np.array, float or double
the real part of the matrix vector multiplication per spin
vcim: 2D np.array, float or double
the imaginary part of the matrix vector multiplication per spin
"""
vcre = np.zeros((2, nspin, nprod), dtype=dtype)
daux = np.zeros((nprod), dtype=dtype)
s2dn = dn.reshape((nspin, nprod))
for spin1 in range(nspin):
for spin2 in range(nspin):
for ireim, sreim in enumerate(('real', 'imag')):
daux[:] = np.require(getattr(s2dn[spin2], sreim), dtype=dtype,
requirements=["A","O"])
vcre[ireim, spin1] += ss2kernel[spin1][spin2].dot(daux) + \
ss2kernel[spin1][spin2].T.dot(daux) - \
ss2kernel_diag[spin1][spin2]*daux
return vcre[0].reshape(-1), vcre[1].reshape(-1)