#* **************************************************************************
#*
#* 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 sys
import math
import numpy as np
import scipy.integrate as spi
import scipy.optimize as spo
[docs]
def modeCouplings(J,rho,n,coord='mf',plot=False,**opts):
"""Return the mode couplings for a bath characterized by its spectral density
This routine implements the equations in Section III of
Wang et al., J. Chem. Phys. 115, 2979 (2001)
J -- function, spectral density of the bath
rho -- discrete bath modes density (Eq. 3.6b)
n -- number of discrete bath modes
coord -- optional, string:
'm': mass-weighted bath coordinates
'mf': mass- and frequency-weighted bath coordinates [DEFAULT]
Returns a list
[ [w1, c1], [w2, c2],....[wn, cn] ]
"""
l = []
fac = 2.0/math.pi
wlist = listOfFrequenciesFromDistribution(n,rho)
for wj in wlist:
cj = math.sqrt( fac*wj*J(wj)/rho(wj) )
if coord=='m':
l.append([wj,cj])
elif coord=='mf':
l.append([wj,cj/math.sqrt(wj)])
# Plot function if requested
if plot:
try:
import matplotlib.pyplot as plt
except ImportError:
message = 'matplotlib module missing. plotting is not possible.\n'
sys.stdout.write(message)
ll = np.array(l)
plt.plot(ll[:,0], ll[:,1],'ro',ll[:,0], ll[:,1],'b-')
plt.show()
return l
[docs]
def listOfFrequenciesFromDistribution(n,rho):
"""Return list of freqs at which Eq 3.6b is satisfied for n modes
This condition means \int_0^w rho(w)dw = k
n -- number of bath modes
rho -- number of bath modes within a frequency interval (Eq. 3.6b)
"""
l=[]
for j in range(n):
k = j + 1
wk = _rho_min(k,rho)
l.append(wk)
return l
def _rho_min(k,rho):
def _rho_error(w):
return (spi.quad(rho,0,w)[0] - k )**2
result = spo.minimize(_rho_error,0)
return result.x[0]
[docs]
def spectralDensityOhmic(alpha,wc,cutoff='sharp',plot=False,**opts):
"""Return an Ohmic spectral density
J(w) = (pi/2) * alpha * w * cutoff(w)
alpha -- ohmic coupling constant
wc -- cutoff frequency
cutoff -- string
'sharp': cutoff(w) = Heaviside(w-wc)
'exponential': cutoff(w) = exp(-w/wc)
plot -- logical, plot the function if set to true
"""
pi2 = math.pi/2.0
# Define spectral density function
if cutoff == 'sharp':
def J(w):
if w > wc: return 0.0
return pi2*alpha*w
elif cutoff == 'exponential':
def J(w):
return pi2*alpha*w*math.exp(-w/wc)
else:
raise ValueError('cutoff == '+cutoff+' unknown')
# Plot function if requested
if plot:
try:
import matplotlib.pyplot as plt
except ImportError:
message = 'matplotlib module missing. plotting is not possible.\n'
sys.stdout.write(message)
npoints = opts.get('xpoints',500)
xfactor = opts.get('xfactor',5)
x = np.linspace(0,wc*xfactor,npoints)
y = np.array(map(J,x))
plt.plot(x, y, 'g-', linewidth=2.0)
plt.xlabel('omega')
plt.ylabel('J')
plt.show()
return J
[docs]
def bathModeDensity(N,cut,shape='constant',wc=1.0):
"""Return a discrete bath mode density function
N -- number of discrete modes in bath
cut -- largest frequency mode
shape -- string, shape of modes distribution
'constant': rho(w) = N/cut
'exponential': rho(w) = (N/wc) * exp(-w/wc) / (1 - exp(-cut/wc))
The function returned is such that the condition
int_0^cut rho(w)dw == N
is always met.
"""
if shape == 'constant':
def rho(w):
return N/cut
elif shape == 'exponential':
# calculate normalization constant
def uexp(x):
return math.exp(-x)
c = N/spi.quad(uexp,0,cut)[0]
print(c)
def rho(w):
return c*math.exp(-w)
else:
raise ValueError('shape == '+shape+' unknown')
return rho