#* **************************************************************************
#*
#* 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.Models.DiabaticToAdiabatic2States import DiabaticToAdiabatic2States as di2ad
# ----------------------------------------------------------------------------
# Tully's model A (single avoided crossing)
# JCP 93, 1061 (1990)
# ----------------------------------------------------------------------------
# Diabatic model
#@np.vectorize
[docs]
def W11_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
if x >= 0:
return A*(1.0-np.exp(-B*x))
if x < 0:
return -A*(1.0-np.exp(B*x))
#@np.vectorize
[docs]
def W22_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return -W11_TA(x,A,B,C,D)
#@np.vectorize
[docs]
def W12_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return C*np.exp(-D*x**2)
#@np.vectorize
[docs]
def W_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
W = np.empty((2,2),float)
W[0,0] = W11_TA(x,A,B,C,D)
W[1,1] = -W[0,0]
W[0,1] = W12_TA(x,A,B,C,D)
W[1,0] = W[0,1]
return W
#@np.vectorize
[docs]
def dq_W11_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
# Derivative of diabatic matrix element wrt all coordinates (only one in this model)
# All derivatives are returned by convention in a list
if x >= 0:
return np.array( (A*B*np.exp(-B*x) ),float )
if x < 0:
return np.array( (A*B*np.exp(B*x) ),float )
#@np.vectorize
[docs]
def dq_W22_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return -1.0*dq_W11_TA(x,A,B,C,D)
#@np.vectorize
[docs]
def dq_W12_TA(x,A=0.01,B=1.6,C=0.005,D=1.0):
return np.array( ( -C*np.exp(-D*x**2)*2.0*D*x ),float )
mod=di2ad(W11_TA,W22_TA,W12_TA,dq_W11_TA,dq_W22_TA,dq_W12_TA)
v1 = mod.ret_V1()
v2 = mod.ret_V2()
dqv1 = mod.ret_dqV1()
dqv2 = mod.ret_dqV2()
d12 = mod.ret_NAC()
# ----------------------------------------------------------------------------
# Interface arguments. Names matter!
# ----------------------------------------------------------------------------
[docs]
def f_EGrad(x,v,t,s):
if s==0:
return v1(x),dqv1(x)
elif s==1:
return v2(x),dqv2(x)
else:
raise ValueError
[docs]
def f_NAC(x,v,t,i,j):
return d12(x)
[docs]
def f_W(x,v,t):
W = np.zeros((2,2),float)
W[0,0] = v1(x)
W[1,1] = v2(x)
return W
[docs]
def f_D(x,v,t):
D = np.zeros((2,2),complex)
D[0,1] = -1.0j*np.dot( d12(x), v )
D[1,0] = np.conj(D[0,1])
return D
masses = [2000.0,]