#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* 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/>.
#*
#* **************************************************************************
#!/usr/bin/python
import sys
import CDTK.Tools.Conversion as conv
from CDTK.Tools.Mathematics import ovrMatrix2
from math import sin, cos, pi, exp, sqrt
import numpy as np
from numpy import dot
[docs]
class FieldCouplings(object):
"""
General properties of the pulses
General functions provided for the field couplings
"""
def __init__(self,**opts):
self.waveformType = opts.get('waveformType','gaussian') # waveform default [gaussian]. optional [sinuoidal] and [user_defined]
self.waveformKernel = opts.get('waveformKernel',None) # user defined waveform function
self.waveformIntensity = opts.get('waveformIntensity',None) # electric field intensity
self.unitIntensity = opts.get('unitIntensity','SI') # unit of electric field intensity. default [SI] and [W_Mm2]
self.waveformAmplitude = opts.get('waveformAmplitude',None) # electric field amplitude
self.unitAmplitude = opts.get('unitAmplitude','atomic_unit') # unit of electric field amplitude. default [atomic_unit]
self.waveformFrequency = opts.get('waveformFrequency',None) # wave frequency
self.unitFrequency = opts.get('unitFrequency','atomic_unit') # unit of wave frequency
self.polarisation = opts.get('polarisation',None) # polarisation vector of the electric field
[docs]
def getElectricField(self,time,**opts):
"""
Return the current electric field strength
"""
if self.waveformType == 'gaussian':
efield = self.getGaussianPulse(time)
elif self.waveformType == 'sinusoidal':
efield = self.getSinusoidalPulse(time)
elif self.waveformType == 'user_defined':
efield = self.getUserDefinedPulse(time)
else:
raise ValueError('waveform of electric field not defined')
return efield
[docs]
def getGaussianPulse(self,time,**opts):
"""
Gaussian pulse
"""
self.convert_waveform_unit()
pass
[docs]
def getSinusoidalPulse(self,time,**opts):
"""
Sinusoidal pulse
"""
self.convert_waveform_unit()
amplitude = self.waveformAmplitude
frequency = self.waveformFrequency
polarisation = self.polarisation
pulse = amplitude * sin(frequency * time) * np.array(polarisation)
return pulse
[docs]
def getUserDefinedPulse(self,time,**opts):
"""
User defined pulse
Pulse [Px Py Pz][TIME]
"""
polarisation = self.polarisation
pulse = self.waveformKernel(time) * np.array(polarisation)
return pulse
[docs]
def getEpF_D_Umat(self, x, nstates, dxgrad, external_E, MuE_MatElt):
"""
Return the gradient of the transformation matrix [EpF_D_Umat] of adiabatic-dressed representations
"""
dx = dxgrad # step for numerical differentiation
EpF_D_Umat = np.zeros((x.size, nstates, nstates),float)
for i in range(x.size):
x1 = x.copy()
x1[i] = x1[i] + 0.5*dx
e1_adiabatic = np.diag(external_E(x1)) - MuE_MatElt
e1_dressed, Umat1 = np.linalg.eig(e1_adiabatic)
x0 = x.copy()
x0[i] = x0[i] - 0.5*dx
e0_adiabatic = np.diag(external_E(x0)) - MuE_MatElt
e0_dressed, Umat0 = np.linalg.eig(e0_adiabatic)
EpF_D_Umat[i,:,:] = (Umat1 - Umat0) / dx
return EpF_D_Umat
[docs]
def getDressedVd(self, Vd, v0, nstates, EpF_Umat, EpF_D_Umat):
"""
Returns the dressed non-adiabatic coupling matrix [Vd_dressed]
Implements Eq.(12) in Richter et al. ; JCTC 7, 1253 (2011)
"""
Vd_adiabatic = Vd.copy()
# [Vd_dressed] contains the basis rotation [Vd_basis_rot] and the rotation gradient [Vd_gradient_rot]
# basis rotation
Vd_basis_rot = np.dot(np.transpose(EpF_Umat), np.dot(Vd_adiabatic, EpF_Umat))
# rotation gradient
ncoords = v0.size
UTmat_D_Umat = np.zeros((ncoords, nstates, nstates),float)
for i in range(ncoords):
UTmat_D_Umat[i] = np.dot(np.transpose(EpF_Umat), EpF_D_Umat[i])
Vd_gradient_rot = np.zeros((nstates, nstates),float)
for i in range(ncoords):
Vd_gradient_rot += v0[i] * UTmat_D_Umat[i]
Vd_dressed = Vd_basis_rot + Vd_gradient_rot
return Vd_dressed
[docs]
def getDressedVd_dUdt(self, Vd, nstates, EpF_Umat, EpF_dUdt):
"""
Returns dressed non-adiabatic coupling matrix [Vd_dressed]
calculated from direct dUdt
"""
Vd_adiabatic = Vd.copy()
# [Vd_dressed] contains the basis rotation [Vd_basis_rot] and the rotation gradient [Vd_gradient_rot]
# basis rotation
Vd_basis_rot = np.dot(np.transpose(EpF_Umat), np.dot(Vd_adiabatic, EpF_Umat))
# rotation gradient
Vd_gradient_rot = np.dot(np.transpose(EpF_Umat), EpF_dUdt)
Vd_dressed = Vd_basis_rot + Vd_gradient_rot
return Vd_dressed