#* **************************************************************************
#*
#* 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 math
import numpy as np
"""
Simple avoided crossing model potential
"""
[docs]
class SimpleACM(object):
def __init__(self,**opts):
self.natoms = opts.get('natoms',1)
self.ndim = opts.get('ndim',3)
self.nstates = opts.get('nstates',2)
if not self.natoms == 1:
raise AssertionError('Only single atom available currently for the model')
def _potential(self,x,**opts):
"""
Definition of the potential of simple avoided crossing
"""
# potential energy
retAll = opts.get('returnAll',False)
retstats = opts.get('returnStates',[])
# cpl matrix element
retAMat = opts.get('returnAMat',False)
xx = np.asarray(x)
xx.shape = (self.natoms,self.ndim)
null = np.array([0],float)
rij = xx[0] - null
Rij = math.sqrt(np.dot(rij,rij))
# potential parameters
A = 0.01
B = 1.60
C = 0.005
D = 1.00
if Rij >= 0.0:
w11 = A * ( 1.0 - np.exp( -B * Rij ) )
w22 = -w11
elif Rij < 0.0:
w11 = -A * ( 1.0 - np.exp( B * Rij ) )
w22 = -w11
w12 = C * np.exp( -D * Rij * Rij )
w21 = w12
# potential in diabatic representation
wdia = np.array([w11, w12, w21, w22], float)
wdia.shape = (self.nstates,self.nstates)
# diagonalize W for adiabatic representation
eval, evec = np.linalg.eig(wdia)
xx.shape = (self.ndim*self.natoms,)
wdia.shape = (self.nstates*self.nstates,)
if retstats:
Epes = []
for s in retstats:
pes = eval[s]
Epes.append(pes)
return np.array(Epes)
elif retAll:
Epes = []
# currently restricted to two electronic potentials
idx = np.arange(self.nstates, dtype=int)
for i in idx:
pes = eval[i]
Epes.append(pes)
return np.array(Epes), wdia
if retAMat:
return evec
def _cplmatrixelement(self,x):
"""
Cpl matrix element d12
"""
dx = 0.0001
Mat_A = np.linalg.inv( self._potential( x, returnAll=False, returnAMat=True ) )
x_v = x + np.array( [dx], float )
Mat_A_v = np.linalg.inv( self._potential( x_v, returnAll=False, returnAMat=True ) )
x_r = x - np.array( [dx], float )
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)
cpl_D = abs(-np.dot(dA,A_inv))
cpl_D.shape = (self.nstates*self.nstates, self.natoms*self.ndim)
for i, cpl_D_matelt in enumerate(cpl_D):
cpl_D[i] = np.array([cpl_D_matelt],float)
return cpl_D
[docs]
def getPotentialFunction(self):
"""
Interface of potential
"""
def f(x):
return self._potential(x, returnAll=True, returnCME=False)
return f
[docs]
def getCplMatrixElement(self):
"""
Interface of coupling matrix elements
"""
def c(x):
return self._cplmatrixelement(x)
return c