# * **************************************************************************
# *
# * CDTK, Chemical Dynamics Toolkit
# * A modular system for chemical dynamics applications and more
# *
# * Copyright (C) 2011, 2012, 2013, 2014, 2015
# * Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
# * Ralph Welsch, DESY, <ralph.welsch@desy.de>
# *
# * Copyright (C) 2020, 2021, 2022, 2023
# * Ludger Inhester, DESY, ludger.inhester@cfel.de
# *
# * This file is part of CDTK.
# *
# * CDTK is free software: you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation, either version 3 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program. If not, see <http://www.gnu.org/licenses/>.
# *
# * **************************************************************************
import sys
import numpy as np
import CDTK.Tools.Utils as uti
from CDTK.Tools.Utils import getUniqueID
import CDTK.Tools.EField as efi
import CDTK.Tools.Conversion as cv
from scipy import interpolate
import numpy.linalg as la
import math
TINY = 1.e-8
[docs]
class interpolation(object):
def __init__(self, **opts):
self.engine = 'interpolation'
self.filepath_energy = opts.get('filepath_energy', None)
self.filepath_td = opts.get('filepath_td', None)
self.filepath_nac = opts.get('filepath_nac', None)
# No. of states + the atomic distance as the a column
self.nums_r = opts.get('nums_r', 2)
self.e_ref = opts.get('e_ref', 0.0)
self.atomicSymbols = opts.get('atomicSymbols', [])
self.atomicNumbers = opts.get('atomicNumbers', [])
self.efields = opts.get('efields', [])
self._X = 1.0
[docs]
def f_E_Grd(self, X, t):
"""
Returns interpolated energy, gradient and Transistion dipole element in a.u. of the current state.
Energy as an one dimesional array and the gradients as a n*6 dimension array.
X -- 3N array with atomic coordinates
S -- electronic state. This is disregarded here
t -- time
"""
a = X[0:3]-X[3:6]
R = la.norm(a)
self._X = R
f = open(str(self.filepath_energy), 'r')
lines = f.readlines()
x = np.linspace(0, 8, self.nums_r)
c = {}
for i in range(len(x)):
c[i] = []
for line in lines:
for j in range(len(x)):
p = line.split()
if j == 0:
c[j].append(float(p[j]))
else:
c[j].append(float(p[j])+float(self.e_ref))
tck = {}
for i in range(len(x)-1):
tck[i] = []
e = []
g = []
for k in range(len(x)-1):
tck[k] = interpolate.splrep(c[0], c[k+1], s=0)
e.append(interpolate.splev(R, tck[k], der=0))
g.append(interpolate.splev(R, tck[k], der=1))
g_N = []
xx = [0, 1]
xxx = [0, 1, 2]
for k in range(len(g)):
for kk in xx:
for kkk in xxx:
if kk > 0:
g_N.append(-(g[k]*a[kkk])/R)
else:
g_N.append((g[k]*a[kkk])/R)
if len(x)-1 > 1:
g_N = np.asarray(g_N).flatten().reshape((len(x)-1, 6))
else:
g_N = np.asarray(g_N).flatten()
e = np.asarray(e)
return e, g_N
[docs]
def f_TD(self, X, t):
a = X[0:3]-X[3:6]
R = la.norm(a)
self._X = R
x = np.linspace(0, 8, self.nums_r)
f = open(str(self.filepath_energy), 'r')
lines = f.readlines()
c = {}
for i in range(len(x)):
c[i] = []
for line in lines:
for j in range(len(x)):
p = line.split()
if j == 0:
c[j].append(float(p[j]))
else:
c[j].append(float(p[j])+float(self.e_ref))
tck = {}
for i in range(len(x)-1):
tck[i] = []
e = []
for k in range(len(x)-1):
tck[k] = interpolate.splrep(c[0], c[k+1], s=0)
e.append(interpolate.splev(R, tck[k], der=0))
e = np.asarray(e)
if self.filepath_td:
teta = np.arccos((R*abs(a[2]))/(R*R))
Rot = np.array([[np.cos(teta), 0.0, np.sin(teta)], [
0.0, 1.0, 0.0], [-np.sin(teta), 0.0, np.cos(teta)]], np.float64)
h = open(str(self.filepath_td), 'r')
lines1 = h.readlines()
dummy = {}
for ii in range(((len(x)-1)*(len(x)-1)*3)+1):
dummy[ii] = []
for l in lines1:
for ii in range(((len(x)-1)*(len(x)-1)*3)+1):
p1 = l.split()
dummy[ii].append(float(p1[ii]))
intr_td = {}
for i in range((len(x)-1)*(len(x)-1)*3):
intr_td[i] = []
td = []
for jj in range((len(x)-1)*(len(x)-1)*3):
intr_td[jj] = interpolate.splrep(dummy[0], dummy[jj+1], s=0)
td.append(interpolate.splev(R, intr_td[jj], der=0))
mu = {}
for ii in range((len(x)-1)):
for jj in range((len(x)-1)):
mu[ii, jj] = []
for ii in range((len(x)-1)):
for jj in range((len(x)-1)):
mu[ii, jj].extend(
[td[(3*jj)+0+(24*ii)], td[(3*jj)+1+(24*ii)], td[(3*jj)+2+(24*ii)]])
mu[ii, jj] = np.dot(Rot, np.asarray(mu[ii, jj])).flatten()
else:
s = (1, 3)
mu = {}
for ii in range(len(x)-1):
for jj in range(len(x)-1):
mu[ii, jj] = np.zeros(s)
# print mu[0,2], R*cv.au2an
return e, mu
[docs]
def f_NAC(self, X, t):
a = X[0:3]-X[3:6]
R = la.norm(a)
x = np.linspace(0, 8, self.nums_r)
f = open(self.filepath_energy, 'r')
lines = f.readlines()
c = {}
for i in range(len(x)):
c[i] = []
for line in lines:
for j in range(len(x)):
p = line.split()
if j == 0:
c[j].append(float(p[j]))
else:
c[j].append(float(p[j])+float(self.e_ref))
tck = {}
for i in range(len(x)-1):
tck[i] = []
e = []
for k in range(len(x)-1):
tck[k] = interpolate.splrep(c[0], c[k+1], s=0)
e.append(interpolate.splev(R, tck[k], der=0))
e = np.asarray(e)
if self.filepath_nac:
# TODO: have to code it.
s = (1, 3)
nac = {}
for ii in range(len(x)-1):
for jj in range(len(x)-1):
nac[ii, jj] = np.zeros(s)
else:
s = (1, 6)
nac = {}
for ii in range(len(x)-1):
for jj in range(len(x)-1):
nac[ii, jj] = np.zeros(s).flatten()
return e, nac
[docs]
def getPartialCharges(self):
return np.zeros(2)
[docs]
def getGamma(self):
"""
Gives the auger rates which is required for the large hops.
"""
probs = [0.00020409336568968513, 0.0014362125733718583, 0.00026456547404218444,
0.00071810628668592916, 0.0004157457449234327, 0.0, 0.0, 0.00035527363657093336]
gamma = []
X = np.zeros(6)
X[5] = self._X
e, g = self.f_E_Grd(X, 0.0)
file_old = self.filepath_energy
num_old = self.nums_r
self.filepath_energy = 'PEC.data'
self.nums_r = 9
e_i, g = self.f_E_Grd(X, 0.0)
self.filepath_energy = file_old
self.nums_r = num_old
for i in range(8):
gamma.append(('a', str(i), probs[i], e[0]-e_i[i]))
return gamma
# [('a','0',0.00020409336568968513,0),('a','1',0.0014362125733718583,0),('a','2',0.00026456547404218444,0),('a','3',0.00071810628668592916,0),('a','4',0.0004157457449234327,0),('a','5',0.0,0),('a','6',0.0,0),('a','7',0.00035527363657093336,0)]
[docs]
def get_occ(self):
return None
[docs]
def set_occ(self, occ):
return None
[docs]
def pe(self, t=0.0):
x = 0.0
return x
[docs]
def efield(self, t=0.0):
"""
Return array with x,y,z components of E field
"""
E = np.zeros(3, float)
for ef in self.efields:
E = E + ef.P*ef.A(t)
return E
[docs]
def dt_efield(self, t=0.0, Dt=0.1):
"""
Return array with x,y,z components of dE/dt (field time derivative)
"""
dEdt = np.zeros(3, float)
for ef in self.efields:
Em = ef.P*ef.A(t-0.5*Dt)
Ep = ef.P*ef.A(t+0.5*Dt)
dEdt = dEdt + (Ep-Em)/Dt
return dEdt