#* **************************************************************************
#*
#* 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 as mt
import numpy as np
[docs]
class DiabaticToAdiabatic2States(object):
"""
Provide adiabatic quantities from a 2-states diabatic model Hamiltonian
"""
def __init__(self,fW11,fW22,fW12,fdqW11,fdqW22,fdqW12):
self.fW11 = fW11
self.fW22 = fW22
self.fW12 = fW12
self.fdqW11 = fdqW11
self.fdqW22 = fdqW22
self.fdqW12 = fdqW12
[docs]
def ret_V1(self):
"""
Return pot function for lower eigenstate
"""
def f(X):
w11 = self.fW11(X)
w22 = self.fW22(X)
w12 = self.fW12(X)
E = w11 + w22
D = w11 - w22
return 0.5*( E - np.sqrt( D**2 + 4*w12**2) )
return f
[docs]
def ret_V2(self):
"""
Return pot function for upper eigenstate
"""
def f(X):
w11 = self.fW11(X)
w22 = self.fW22(X)
w12 = self.fW12(X)
E = w11 + w22
D = w11 - w22
return 0.5*( E + np.sqrt( D**2 + 4.0*w12**2) )
return f
[docs]
def ret_NAC(self):
"""
Return function for non-adiabatic coupling vector
"""
def f(X):
w11 = self.fW11(X)
w22 = self.fW22(X)
w12 = self.fW12(X)
dq_w11 = self.fdqW11(X)
dq_w22 = self.fdqW22(X)
dq_w12 = self.fdqW12(X)
Q = 2.0*w12/( w11 - w22 )
P = 0.5*np.arctan(Q)
DE = w11 - w22
dq_DE = dq_w11 - dq_w22
dq_P = 1.0/(1.0+Q**2)*( dq_w12*DE - w12*dq_DE )/DE**2
return dq_P
return f
[docs]
def ret_dqV1(self):
"""
Gradient of lower eigenstate
"""
def f(X):
w11 = self.fW11(X)
w22 = self.fW22(X)
w12 = self.fW12(X)
dq_w11 = self.fdqW11(X)
dq_w22 = self.fdqW22(X)
dq_w12 = self.fdqW12(X)
E = w11 + w22
D = w11 - w22
R = np.sqrt( D**2 + 4*w12**2)
dqE = dq_w11 + dq_w22
dqD = dq_w11 - dq_w22
return 0.5*( dqE - 0.5/R*( 2.0*D*dqD + 8.0*w12*dq_w12 ) )
return f
[docs]
def ret_dqV2(self):
"""
Gradient of upper eigenstate
"""
def f(X):
w11 = self.fW11(X)
w22 = self.fW22(X)
w12 = self.fW12(X)
dq_w11 = self.fdqW11(X)
dq_w22 = self.fdqW22(X)
dq_w12 = self.fdqW12(X)
E = w11 + w22
D = w11 - w22
R = np.sqrt( D**2 + 4*w12**2)
dqE = dq_w11 + dq_w22
dqD = dq_w11 - dq_w22
return 0.5*( dqE + 0.5/R*( 2.0*D*dqD + 8.0*w12*dq_w12 ) )
return f