#* **************************************************************************
#*
#* 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
"""
Harmonic oscillator model potential
"""
[docs]
class HarmonicOsc(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',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 harmonic oscillators
Options -
enable_dipole : whether to calculate dipole moment
"""
# potential energy
#retAll = opts.get('returnAll',False)
retAll = opts.get('returnAll',True)
retstats = opts.get('returnStates',[])
#enable_dipole = opts.get('enable_dipole',False)
enable_dipole = opts.get('enable_dipole',True)
xx = np.asarray(x)
xx.shape = (self.natoms,self.ndim)
null = np.zeros(self.ndim*self.natoms,dtype=float)
rij = xx[0] - null
Rij = math.sqrt(np.dot(rij,rij))
# potential parameters
k = 1.0
D12 = 40.0
v11 = 0.5 * k * Rij * Rij
v22 = 0.5 * k * Rij * Rij + D12
# potential in adiabatic representation
vadia = np.array([v11, v22])
xx.shape = (self.ndim*self.natoms,)
# dipole moment matrix
# each matrix element is one-dimensional vector
dipole = []
for i in range(self.nstates):
dipole.append([])
for j in range(self.nstates):
if j != i:
dipole_x = [1.0]
else:
dipole_x = [0.0]
dipole[i].append(dipole_x)
if retstats:
Epes = []
for s in retstats:
pes = vadia[s]
Epes.append(pes)
if enable_dipole is True:
return np.array(Epes), np.array(dipole)
else:
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 = vadia[i]
Epes.append(pes)
if enable_dipole is True:
return np.array(Epes), np.array(dipole)
else:
return np.array(Epes)
def _cplmatrixelement(self,x):
"""
Cpl matrix element
"""
cpl_D = np.zeros([self.nstates,self.nstates],float)
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,**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):
return self._cplmatrixelement(x)
return c