#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
import numpy as np
from CDTK.Tools.Mathematics import sortEigValsVecs as sort
"""
Multiple diabats model potential
"""
[docs]
class MultDiabats(object):
def __init__(self,**opts):
"""
natoms : number of atoms
ndim : spatial dimensions
nstates : number of electronic states
"""
self.natoms = opts.get('natoms',1)
self.ndim = opts.get('ndim',1)
self.nstates = opts.get('nstates',2)
if not self.natoms == 1:
raise AssertionError('Only single atom available currently for the model')
if not self.ndim == 1:
raise AssertionError('Only one dimension available currently for the model')
def _potential(self,x,**opts):
"""
Definition of the potential of harmonic oscillators
Options -
enable_dipole : whether to calculate dipole moment
"""
# potential energy
# adiabat
retAll = opts.get('returnAll',False)
retstats = opts.get('returnStates',[])
# diabat
retDia = opts.get('returnDiabats',False)
# cpl matrix element
retAMat = opts.get('returnAMat',False)
xx = x.copy()
xx.shape = (self.natoms,self.ndim)
null = np.zeros(self.ndim*self.natoms,dtype=float)
X = xx[0,0]
# potential parameters
kappa = 0.50 # slope of the crossing curve
delta = 0.50 # energy spacing of the parallel curves
gamma = 0.05 # off-diagonal coupling of diabats
wdia = np.zeros([self.nstates, self.nstates], dtype=float)
# state zero
wdia[0,0] = - kappa * X
# state one to the roof
E_roof = 0.0
for i in range(1, self.nstates):
E_s = (i - 1) * delta # energy origin shift downwards
# diagonal elements
wdia[i,i] = E_roof - E_s
# off-diagonal elements
wdia[i,0] = gamma
wdia[0,i] = wdia[i,0]
# diagonalize W for adiabatic representation
eval, evec = np.linalg.eig(wdia)
sort(eval, evec)
xx.shape = (self.ndim*self.natoms,)
if retstats:
Epes = []
for s in retstats:
pes = eval[s]
Epes.append(pes)
return np.array(Epes)
elif retAll:
Epes = []
idx = np.arange(self.nstates, dtype=int)
for i in idx:
pes = eval[i]
Epes.append(pes)
return np.array(Epes)
elif retDia:
Epes = []
idx = np.arange(self.nstates, dtype=int)
for i in idx:
pes = wdia.diagonal()[i]
Epes.append(pes)
return np.array(Epes)
if retAMat:
return evec
def _cplmatrixelement(self,x,state_block):
"""
Non-adiabatic coupling matrix element
"""
dx = 0.001
D = []
for i in range(self.ndim):
Mat_A = np.linalg.inv( self._potential( x, returnAll=False, returnAMat=True ) )
xx = x.copy()
xx[i] += dx
x_v = xx.copy()
Mat_A_v = np.linalg.inv( self._potential( x_v, returnAll=False, returnAMat=True ) )
xx = x.copy()
xx[i] -= dx
x_r = xx.copy()
Mat_A_r = np.linalg.inv( self._potential( x_r, returnAll=False, returnAMat=True ) )
dA = ( ( Mat_A_v - Mat_A_r ) / ( 2.0 * dx ) )
A_inv = np.linalg.inv(Mat_A)
D.append( - np.dot(dA,A_inv) )
# construct the blockwise coupling matrix
for D_elt in D:
for i in range(self.nstates):
for j in range(self.nstates):
if ( (not i in state_block) or (not j in state_block) ):
D_elt[i,j] = 0.0
cpl_D = np.zeros((self.nstates*self.nstates, self.natoms*self.ndim),dtype=float)
for i, D_elt in enumerate(D): # idx for natoms x ndim
for j in range(self.nstates):
for k in range(self.nstates):
cpl_D[j*self.nstates+k,i] = D_elt[j,k]
for i, cpl_D_matelt in enumerate(cpl_D):
cpl_D[i] = np.array([abs(cpl_D_matelt)],float)
return cpl_D
[docs]
def getPotentialFunction(self,**opts):
"""
Interface of potential
"""
def f(x,**opts):
return self._potential(x,**opts)
return f
[docs]
def getCplMatrixElement(self):
"""
Interface of coupling matrix elements
"""
def c(x, state_block):
return self._cplmatrixelement(x,state_block)
return c