from __future__ import print_function, division
import numpy as np
from timeit import default_timer as timer
from pynao.m_siesta_units import siesta_conv_coefficients
"""
Iterative TDDFT using the electostatic potential of a moving charge as
perturbation. The units of the input are in atomic units.
The theory behind this implementation and some application can be found
in chapter 6 of [Barbry2018thesis]_.
"""
[docs]def check_collision(velec, beam_offset, atom2coord, threshold=1.0e-6, verbosity=1):
"""
Check if the electron collide with an atom
if the trajectory of the electron is passing too close of an atom, i.e.,
if at any moment in time the distance between the electron and an atom
is below the value given by threshold.
Also check that the electron velocity and the beam ofset are perpendicular.
Input Parameters:
atom2coord: 2D np.array, float
Atoms coordinates
threshold: float
Threshold use to determine if the beam is colliding an atom.
"""
vnorm = np.sqrt(np.dot(velec, velec))
vdir = velec/vnorm
if verbosity>0:
print("{} => tem parameters:".format(__name__))
print("vdir: ", vdir)
print("vnorm: ", vnorm)
print("beam_offset: ", beam_offset)
# check orthogonality between beam direction
# and beam offset
assert abs(np.dot(velec, beam_offset)) < threshold
R0 = -100.0*np.max(atom2coord)*vdir + beam_offset
for atm in range(atom2coord.shape[0]):
vec = R0 - atom2coord[atm, :]
# unit vector to compare to vdir
vec = abs(vec/np.sqrt(np.dot(vec, vec)))
if np.sqrt(np.dot(vec - vdir, vec - vdir)) < threshold:
######### fancy message does not work in python2
mess = 'np.sqrt(np.dot(vec - vdir, vec - vdir))<1e-6:'
print("atoms {0} coordinate: ".format(atm), atom2coord[atm, :])
#mess = """
#Electron is collinding with atom {0}:
#velec = [{1:.3f}, {2:.3f}, {3:.3f}]
#beam_offset = [{4:.3f}, {5:.3f}, {6:.3f}]
#atom coord = [{7:.3f}, {8:.3f}, {9:.3f}]
#impact parameter = {10:.9f} > 1e-6""".format(atm, *self.velec,
# *self.beam_offset[0],
# *atom2coord[atm, :],
# np.sqrt(np.dot(vec, self.vdir)))
raise ValueError(mess)
[docs]def get_time_range(velec, beam_offset, dr, freq, verbosity=0):
"""
Get the time and symmetric frequency range for the electron passing close
to the particle. The time range is a symmetric array
around 0.0. At t = 0, the electron is at its closest
position from the molecule. This array will depend on the
frequency range and the spatial precision dr.
To respect the Fourier transform convention, the following
relationshiip must be fulfill,
N = 2*pi/(dw*dt)
with N the number of element of t.
N must be an odd number in order that t is symmetric
Input Parameters:
dr: 1D np.array, float, size 3
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 highly inacurate.
freq: 1D np.array, float
Frequency range (in atomic unit), freq[0] must be 0.0!!
"""
from pynao.m_tools import is_power2
vnorm = np.sqrt(np.dot(velec, velec))
vdir = velec/vnorm
dt = np.min(dr)/vnorm
dw = freq[1] - freq[0]
N_org = int(2*np.pi/(dw*dt))
# to improve performance, N must be a power of 2
if not is_power2(N_org):
power = 1
while 2**power < N_org:
power +=1
minima = np.argmin(np.array([abs(2**(power-1) - N_org), abs(2**power - N_org)]))
if minima == 0:
N = 2**(power-1)
else:
N = 2**power
if verbosity > 0:
print("N_org = {0}, N_new = {1}".format(N_org, N))
dt = 2*np.pi/(N*dw)
dr = np.array([dt*vnorm, dt*vnorm, dt*vnorm])
else:
N = N_org
dw_symm = 2.0*np.pi/(N*dt)
wmax = 2.0*np.pi*(N-1)/(N*dt)/2.0
freq_symm = np.arange(-wmax, wmax+dw_symm, dw_symm)[0:N]
tmax = (N-1)*dt/2
time_range = np.arange(-tmax, tmax+dt, dt)[0:N]
return time_range, freq_symm
[docs]def calc_external_potential(nprod, velec, beam_offset, freq, freq_sym, time_range,
ao_log=None, verbosity=0):
"""
Calculate the external potential created by a moving charge
"""
from pynao.m_comp_vext_tem import comp_vext_tem
t1 = timer()
V_freq = comp_vext_tem(nprod, velec, beam_offset, freq, freq_sym, time_range,
ao_log=ao_log)
t2 = timer()
if verbosity > 2:
print("{} => calc_external_potential timing: {}".format(__name__, t2-t1))
print("sum(V_freq) = ", np.sum(abs(V_freq.real)),
np.sum(abs(V_freq.imag)))
return V_freq
[docs]def comp_tem_spectrum(self, comegas, V_freq, tmp_fname=None, inter=True):
"""
Compute the interacting tem spectrum
Input Parameters:
comegas: 1D np.array, complex
frequency range (in Hartree) for which the polarizability is computed.
The imaginary part control the width of the signal.
For example,
>>> td = tddft_iter_c(...)
>>> comegas = np.arange(0.0, 10.05, 0.05) + 1j*td.eps
V_freq: 2D np.array, complex
The perturbation to the system
inter: boolean
Perform TDDFT calculation without or with interaction
tmp_fname: string
temporary file to save polarizability at each frequency.
Can be a life saver for large systems.
The format of the file is the following,
# energy (Hartree) Re(gamma) Im(gamma)
Output Parameters:
dn: 2D np.array, complex
computed density change in prod basis
gamma: 1D np.array, complex
computed eels spectrum
"""
gamma = np.zeros_like(comegas, dtype=self.dtypeComplex)
dn = np.zeros((comegas.shape[0], self.nprod), dtype=self.dtypeComplex)
for iw, comega in enumerate(comegas):
chi0mv_ncalls_ite = self.rf0_ncalls
t1 = timer()
if inter:
veff = self.comp_veff(V_freq[iw, :], comega)
dn[iw, :] = self.apply_rf0(veff, comega, self.chi0_mv)
gamma[iw] = np.dot(np.conj(V_freq[iw, :]), dn[iw, :])
else:
dn[iw, :] = self.apply_rf0(V_freq[iw, :], comega, self.chi0_mv)
gamma[iw] = np.dot(dn[iw, :], np.conj(V_freq[iw, :]))
t2 = timer()
chi0mv_ncalls_ite = self.rf0_ncalls - chi0mv_ncalls_ite
if self.verbosity > 1:
mess = "iw: {0}/{1}; w: {2:.4f};".format(iw, comegas.shape[0],
comega.real*siesta_conv_coefficients["ha2ev"])
mess += " nite: {0}; time: {1}".format(chi0mv_ncalls_ite, t2-t1)
print(mess)
if tmp_fname is not None:
tmp = open(tmp_fname, "a")
tmp.write("{0} {1} {2}\n".format(comega.real, -gamma[iw].real/np.pi,
-gamma[iw].imag/np.pi))
tmp.close()
# Need to open and close the file at every freq, otherwise
# tmp is written only at the end of the calculations, therefore,
# it is useless
if inter:
self.write_chi0_mv_timing("tddft_tem_inter_chi0_mv.txt")
else:
self.write_chi0_mv_timing("tddft_tem_nonin_chi0_mv.txt")
print("Total number of iterations: ", self.rf0_ncalls)
return dn, -gamma/np.pi