from __future__ import print_function, division
from copy import copy
import numpy as np
import scipy.sparse as sparse
from pynao.chi0_matvec import chi0_matvec
from pynao.m_chi0_noxv import chi0_mv_gpu, chi0_mv
import pynao.m_kernel_utils as kutils
[docs]class tddft_iter(chi0_matvec):
"""tddft_iter object.
Iterative TDDFT a la PK, DF, OC JCTC
References:
1. [Barbry2018thesis]_
2. [koval2018pynao]_
3. [koval2010-mbptlcao]_
"""
def __init__(self, **kw):
self.kernel_format = kw['kernel_format'] if 'kernel_format' in kw else "pack"
chi0_matvec.__init__(self, **kw)
self.xc_code_mf = copy(self.xc_code)
if "xc_code" in kw.keys():
self.xc_code = kw['xc_code']
else:
kw['xc_code'] = self.xc_code
if self.GPU:
self.chi0_mv = chi0_mv_gpu
else:
self.chi0_mv = chi0_mv
if not hasattr(self, 'pb'):
print(__name__, 'no pb?')
print(__name__, kw.keys())
raise ValueError("product basis not initialized?")
self.kernel_dim, self.kernel, self.kernel_diag, self.ss2kernel, \
self.ss2kernel_diag = kutils.kernel_initialization(self, **kw)
if self.verbosity > 0:
print(__name__,'\t====> self.xc_code:', self.xc_code)
[docs] def comp_veff(self, vext, comega=1j*0.0, prev_sol=None):
"""
This computes an effective field (scalar potential) given the external
scalar potential.
Solves
Ax = b
with
* b the external potential delta V_ext
* x the effective potential delta V_eff
* A = (1 - K_{Hxc}Chi_0)
"""
cmplx_types = [complex, np.complex64, np.complex128]
iscomplex = False
for cmplxtype in cmplx_types:
if isinstance(comega, cmplxtype):
iscomplex = True
break
if not iscomplex:
raise ValueError("comega of type {} must be complex!".format(comega))
self.matvec_ncalls = 0
nsp = self.nspin*self.nprod
assert len(vext) == nsp, "{} {}".format(len(vext), nsp)
self.comega_current = comega
#if prev_sol is not None:
# preconditionner = self.preconditioning(comega, vext, prev_sol)
#else:
# preconditionner = None
veff_op = sparse.linalg.LinearOperator((nsp, nsp),
matvec=self.vext2veff_matvec,
dtype=self.dtypeComplex)
# right hand side of the equation
rhs = np.require(vext, dtype=self.dtypeComplex, requirements='C')
#print("shape: ", veff_op.shape, rhs.shape)
# Solves Ax = b
#print(self.krylov_solver)
#print("solver options: ", self.krylov_options)
resgm, info = self.krylov_solver(veff_op, rhs, x0=prev_sol,
M=None,
**self.krylov_options)
if info != 0:
print("LGMRES Warning: info = {0}".format(info))
return resgm
[docs] def preconditioning(self, comega, vin, x0):
nsp = self.nspin*self.nprod
rows = np.arange(nsp)
cols = np.arange(nsp)
data = (x0 - vin)/x0
#data = np.zeros((nsp), dtype=self.dtypeComplex)
#data.fill(complex(1.0, 0.0))
#data -= get_diagonal_triangle_matrix(self.kernel, self.kernel_dim)/comega
return sparse.coo_matrix((1.0/data, (rows, cols)),shape=(nsp, nsp),
dtype=self.dtypeComplex)
[docs] def vext2veff_matvec(self, vin):
dn0 = self.apply_rf0(vin, self.comega_current, self.chi0_mv)
vcre, vcim = self.apply_kernel(dn0)
return vin - (vcre + 1.0j*vcim)
[docs] def vext2veff_matvec2(self, vin):
dn0 = self.apply_rf0(vin, self.comega_current, self.chi0_mv)
vcre,vcim = self.apply_kernel(dn0)
return 1 - (vin - (vcre + 1.0j*vcim))
[docs] def apply_kernel(self, dn):
"""
Apply the kernel to chi0_mv, i.e, matrix vector multiplication.
The kernel is stored in pack or pack sparse format.
Input parameters:
dn: 1D np.array, complex
the results of chi0_mv product, i.e, chi0*delta_V_ext
Output parameters:
vcre: 1D np.array, float
The real part of the resulting matvec product kernel.dot(chi0_mv)
vcim: 1D np.array, float
The imaginary part of the resulting matvec product kernel.dot(chi0_mv)
"""
if self.nspin == 1:
if self.kernel_format == "pack":
return kutils.apply_kernel_nspin1_pack(self.spmv, self.nprod,
self.kernel,
dn, dtype=self.dtype)
elif self.kernel_format == "sparse":
return kutils.apply_kernel_nspin1_sparse(self.nprod, self.kernel,
self.kernel_diag,
dn, dtype=self.dtype)
else:
raise ValueError("unknown sparse format {}".format(self.kernel_format))
elif self.nspin == 2:
if self.kernel_format == "pack":
return kutils.apply_kernel_nspin2_pack(self.spmv, self.nprod,
self.nspin, self.ss2kernel,
dn, dtype=self.dtype)
elif self.kernel_format == "sparse":
return kutils.apply_kernel_nspin2_sparse(self.nprod, self.nspin,
self.ss2kernel,
self.ss2kernel_diag,
dn, dtype=self.dtype)
else:
raise ValueError("unknown sparse format {}".format(self.kernel_format))
[docs] def comp_polariz_nonin_xx(self, comegas, tmp_fname=None):
"""
Compute the non-interacting polarizability along the xx direction
"""
self.dn0, self.p0_mat = \
self.comp_dens_along_eext(comegas,
eext=np.array([1.0, 0.0, 0.0]),
tmp_fname=tmp_fname,
inter=False)
return self.p0_mat[0, 0, :]
[docs] def comp_polariz_inter_xx(self, comegas, tmp_fname=None):
"""
Compute the interacting polarizability along the xx direction
"""
self.dn, self.p_mat = \
self.comp_dens_along_eext(comegas,
eext=np.array([1.0, 0.0, 0.0]),
tmp_fname=tmp_fname,
inter=True)
return self.p_mat[0, 0, :]
[docs] def comp_polariz_nonin_ave(self, comegas, tmp_fname=None):
"""
Compute average interacting polarizability
"""
self.dn0, self.p0_mat = \
self.comp_dens_along_eext(comegas,
eext=np.array([1.0, 1.0, 1.0]),
tmp_fname=tmp_fname,
inter=False)
Pavg = np.zeros((self.p0_mat.shape[2]), dtype=self.dtypeComplex)
for i in range(3):
Pavg[:] += self.p0_mat[i, i, :]
return Pavg/3
[docs] def comp_polariz_inter_ave(self, comegas, tmp_fname=None):
"""
Compute average interacting polarizability
"""
self.dn, self.p_mat = \
self.comp_dens_along_eext(comegas,
eext=np.array([1.0, 1.0, 1.0]),
tmp_fname=tmp_fname,
inter=True)
Pavg = np.zeros((self.p_mat.shape[2]), dtype=self.dtypeComplex)
for i in range(3):
Pavg[:] += self.p_mat[i, i, :]
return Pavg/3
polariz_inter_ave = comp_polariz_inter_ave
[docs] def comp_polariz_nonin_edir(self, comegas, eext=np.array([1.0, 1.0, 1.0]),
tmp_fname=None):
"""
Compute average interacting polarizability
"""
self.dn0, self.p0_mat = \
self.comp_dens_along_eext(comegas, eext=eext,
tmp_fname=tmp_fname,
inter=False)
return self.p0_mat
[docs] def comp_polariz_inter_edir(self, comegas, eext=np.array([1.0, 1.0, 1.0]),
tmp_fname=None):
"""
Compute average interacting polarizability
"""
self.dn, self.p_mat = \
self.comp_dens_along_eext(comegas, eext=eext,
tmp_fname=tmp_fname,
inter=True)
return self.p_mat
[docs] def get_eels_spectrum(self, freq, velec=np.array([1.0, 0.0, 0.0]),
beam_offset=np.array([0.0, 0.0, 0.0]),
dr=np.array([0.3, 0.3, 0.3]),
tmp_fname=None, Vext=None, inter=True):
"""
Calculate the interacting TEM spectra for a given electron trajectory.
Typically, can be used to calculate valence Electron Energy Loss
Spectroscopy (eels). The perturbation is created by a moving charge.
See chapter 6 of [Barbry2018thesis]_ for more details.
Input Parameters:
freq: 1D np.array, float
Frequency range for which the eels spectra is calculated in atomic units.
velec: 1D np.array, float
xyz component of the electron velocity in atomic unit
beam_offset: 1D np.array, float
xyz components of the beam offset, must be orthogonal
to velec in atomic unit
dr: 1D np.array, float
spatial resolution for the electron trajectory in atomic unit.
Warning: This parameter influence the accuracy of the calculations.
If taken too large the results will be highly inacurate.
tmp_fname: string
filename to store temporary results of the eels spectra.
if None, the results will not be saved during the run
Vext: 1D np.array, Complex
Optional, if not given, the external potential created by the charge
will be calculated. Can be used to avoid to recalculate Vext if
doing calculation where only the norm of the velocity is changing.
inter: bool
Perform interactive or non-interactive calculations.
Output Parameters:
Vext: 1D np.array, complex
The external potential created by the charge in frequency domain.
Can be used again for further calculations.
eels: 1D np.array, complex
The calculated eels spectrum
The velocity is given in atomic units, it can be converted from energy
(in Hartree) using the relativistic energy kinetic formulation.
Below is a function doing such conversion
Convert electron energy (in eV) to velocity:
>>> def E2vel(E):
>>> # Fine-structure
>>> alpha = 0.0072973525693
>>> num = np.sqrt((alpha**2)*E**2 + 2*E)
>>> den = E*alpha**2 + 1
>>> return num/den
>>> # Convert electron energy in eV to velocity in atomic unit
>>> Ha = 27.211
>>> E = 1e5/Ha # E = 100 keV
>>> v = E2vel(E)
>>> print(v)
"""
from pynao.tddft_tem import get_time_range, comp_tem_spectrum, \
calc_external_potential, check_collision
assert velec.size == 3
assert beam_offset.size == 3
# First element of the frequency range must be 0.0!
assert abs(freq[0]) < 1.0e-12
if tmp_fname is not None:
if not isinstance(tmp_fname, str):
raise ValueError("tmp_fname must be a string")
check_collision(velec, beam_offset, self.atom2coord, verbosity=self.verbosity)
time_range, freq_sym = get_time_range(velec, beam_offset, dr, freq,
verbosity=self.verbosity)
if Vext is None:
# if only the norm of the velocity vector is changing, one can reuse
# the exxternal potential
Vext = calc_external_potential(self.nprod, velec, beam_offset, freq,
freq_sym, time_range,
ao_log=self.pb.prod_log,
verbosity=self.verbosity)
self.dn, eels = comp_tem_spectrum(self, freq + 1.0j*self.eps, Vext,
tmp_fname=tmp_fname, inter=inter)
return Vext, eels