In [1]:
import pickle
import numpy as np
import pycuda.gpuarray as gpuarray
from scipy.special import hyp1f1
import scipy.fftpack as fftpack
import pylab as plt
import scipy.linalg as linalg
#-------------------------------------------------------------------------------------
from pywignercuda_path import SetPyWignerCUDA_Path
SetPyWignerCUDA_Path()
from GPU_Wigner2D_FFT import *
In [2]:
%matplotlib inline
In [3]:
class frame( GPU_FokkerPlank2D_FFT ):
def __init__ (self):
X_gridDIM = 512
P_gridDIM = 512
X_amplitude = 12
P_amplitude = 12
self.mass = 1.
self.omega = 1.
kappa = 1.
dt= 0.005
timeSteps = 20000
skipFrames = 500
mass = 1.
# Diffusion parameter
D_Theta = 0.0143
D_Lambda = 0.000
# Damping parameters (implies another source of diffusion as well)
self.dampingFunction = 'CaldeiraLeggett'
gammaDamping = 0.07
self.epsilon = 0.5
self.grossPitaevskiiCoefficient = 0.0
# Potential and derivative of potential
X2_constant = 0.5*self.mass*self.omega**2
potentialString = '{0}*pow(x,2)'.format(X2_constant)
dPotentialString = '2*{0}*x'.format(X2_constant)
self.fp_Damping_String = ' p*p/sqrt( p*p + {epsilon} ) '.format( epsilon = self.epsilon**2 )
self.SetTimeTrack( dt, timeSteps, skipFrames,
fileName = '/home/rcabrera/DATA/FokkerPlanck2D/X2.hdf5' )
GPU_FokkerPlank2D_FFT.__init__(self,
X_gridDIM,P_gridDIM,X_amplitude,P_amplitude,
kappa,mass,D_Theta,D_Lambda,gammaDamping,potentialString,dPotentialString)
def Set_Initial_Condition_HarmonicOscillator(self,n,omega,x_init,p_init):
"""
Sets self.PsiInitial_XP with the Wigner function of the Harmonic oscillator
"""
self.x_init = x_init
self.p_init = p_init
self.W_init = self.Wigner_HarmonicOscillator(n,omega,x_init,p_init)
In [4]:
instance = frame()
print ' '
print ' Wigner2D propagator with damping '
print ' '
instance.Set_Initial_Condition_HarmonicOscillator (n=0, omega=1 , x_init=0 , p_init=3 )
W = instance.Run( )
In [5]:
instance.LoadEhrenfestFromFile()
In [6]:
print 'Potential'
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot( instance.X_range, instance.Potential(0,instance.X_range) )
ax.set_xlim(-10,10)
ax.set_ylim(-1,60)
ax.set_xlabel('x')
ax.set_ylabel('V')
ax.grid('on')
In [7]:
def PlotWignerFrame( W_input , x_plotRange,p_plotRange):
W = W_input.copy()
W = fftpack.fftshift(W.real)
dp = instance.dP
p_min = -instance.P_amplitude
p_max = instance.P_amplitude - dp
#p_min = -dp*instance.P_gridDIM/2.
#p_max = dp*instance.P_gridDIM/2. - dp
x_min = -instance.X_amplitude
x_max = instance.X_amplitude - instance.dX
global_max = 0.17 # Maximum value used to select the color range
global_min = -0.31 #
print 'min = ', np.min( W ), ' max = ', np.max( W )
print 'final time =', instance.timeRange[-1] ,'a.u. =',\
instance.timeRange[-1]*( 2.418884326505*10.**(-17) ) , ' s '
print 'normalization = ', np.sum( W )*instance.dX*dp
zero_position = abs( global_min) / (abs( global_max) + abs(global_min))
wigner_cdict = {'red' : ((0., 0., 0.),
(zero_position, 1., 1.),
(1., 1., 1.)),
'green' : ((0., 0., 0.),
(zero_position, 1., 1.),
(1., 0., 0.)),
'blue' : ((0., 1., 1.),
(zero_position, 1., 1.),
(1., 0., 0.)) }
wigner_cmap = matplotlib.colors.LinearSegmentedColormap('wigner_colormap', wigner_cdict, 256)
fig, ax = plt.subplots(figsize=(12, 5))
cax = ax.imshow( W ,origin='lower',interpolation='none',\
extent=[ x_min , x_max, p_min, p_max], vmin= global_min, vmax=global_max, cmap=wigner_cmap)
ax.contour(instance.Hamiltonian ,
np.arange(0, 10, 1 ),origin='lower',extent=[x_min,x_max,p_min,p_max],
linewidths=0.25,colors='k')
axis_font = {'size':'24'}
ax.set_xlabel(r'$x$',**axis_font)
ax.set_ylabel(r'$p$',**axis_font)
ax.set_xlim((x_plotRange[0] , x_plotRange[1] ))
ax.set_ylim((p_plotRange[0] , p_plotRange[1] ))
ax.set_aspect(1.)
#ax.grid('on')
cbar = fig.colorbar(cax, ticks=[-0.3, -0.2,-0.1, 0, 0.1, 0.2 , 0.3])
matplotlib.rcParams.update({'font.size': 18})
return fig
In [8]:
plot_init = PlotWignerFrame( instance.W_init.real , (-5.,5) ,(-5,5) )
In [9]:
plot_init = PlotWignerFrame( instance.W_end , (-5.,5) ,(-5,5) )
In [10]:
def PlotMarginals():
W = fftpack.fftshift( instance.W_end.real )
dp = instance.dP
p_min = -instance.P_amplitude
p_max = instance.P_amplitude - dp
W0 = fftpack.fftshift(instance.W_init.real )
marginal_x_init = np.sum( W0 , axis=0 )*dp
marginal_p_init = np.sum( W0 , axis=1 )*instance.dX
marginal_x = np.sum( W, axis=0 )*dp
marginal_p = np.sum( W, axis=1 )*instance.dX
dx = instance.dX
x_min = -instance.X_amplitude
x_max = instance.X_amplitude - instance.dX
#.......................................... Marginal in position
plt.figure(figsize=(10,10))
plt.subplot(211)
plt.plot(instance.X_range, marginal_x_init, '-',label='initial')
plt.plot(instance.X_range, marginal_x, label='final')
#plt.axis([x_min, 0*x_max, -0.01,6])
plt.xlabel('x')
plt.ylabel('Prob')
plt.legend(loc='upper right', shadow=True)
#.......................................... Marginal in momentum
print 'x = ', instance.X_average[0], ' -> ', instance.X_average[-1]
print 'p = ', instance.P_average[0], ' -> ', instance.P_average[-1]
rangeP = np.linspace( p_min, p_max, instance.P_gridDIM )
plt.subplot(212)
plt.plot(rangeP, marginal_p_init ,'-', label='initial')
plt.plot(rangeP, marginal_p , label='final')
plt.axis([p_min, p_max, -0.01, 2])
plt.xlabel('p')
plt.ylabel('Prob')
plt.legend(loc='upper right', shadow=True)
In [11]:
def PlotLogMarginals():
Wend = fftpack.fftshift( instance.W_end ).real
dp = instance.dP
p_min = -instance.P_amplitude
p_max = instance.P_amplitude - dp
W0 = fftpack.fftshift(instance.W_init.real )
marginal_x_init = np.sum( W0 , axis=0 )*dp
marginal_p_init = np.sum( W0 , axis=1 )*instance.dX
marginal_x = np.sum( Wend, axis=0 )*dp
marginal_p = np.sum( Wend, axis=1 )*instance.dX
x_min = -instance.X_amplitude
x_max = instance.X_amplitude - instance.dX
#.......................................... Marginal in position ....................
plt.figure(figsize=(10,10))
plt.subplot(211)
plt.semilogy(instance.X_range, marginal_x_init, '-',label='initial')
plt.semilogy(instance.X_range, marginal_x, label='final')
plt.ylim( 1e-9 , 2)
plt.xlabel('x')
plt.ylabel('Prob')
plt.legend(loc='upper right', shadow=True)
#.......................................... Marginal in momentum
print 'x = ', instance.X_average[0], ' -> ', instance.X_average[-1]
print 'p = ', instance.P_average[0], ' -> ', instance.P_average[-1]
rangeP = np.linspace( p_min, p_max, instance.P_gridDIM )
plt.subplot(212)
plt.semilogy(rangeP, marginal_p_init ,'-', label='initial')
plt.semilogy(rangeP, marginal_p , label='final')
plt.ylim( 1e-9 , 2)
plt.xlabel('p')
plt.ylabel('Prob')
plt.legend(loc='upper right', shadow=True)
In [12]:
PlotMarginals()
In [13]:
PlotLogMarginals()
In [14]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , np.gradient(instance.X_average, instance.dt) , '-',
label = '$\\frac{d}{dt} \\langle x \\rangle $' ,color = 'red', linewidth=1.5)
ax.plot( instance.timeRange , instance.P_average/instance.mass , '--' ,
label='$\\frac{1}{m}\\langle p \\rangle$', linewidth=1.5 )
#ax.set_xlim(0,3.5)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})
ax.set_xlabel('t')
ax.grid();
In [15]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , np.gradient( instance.P_average , instance.dt) ,'-' ,
label = '$\\frac{d}{dt} \\langle p \\rangle $' ,color = 'r' , linewidth=1.5)
ax.plot( instance.timeRange ,
- instance.dPotentialdX_average -2*instance.gammaDamping*instance.P_average , '--' ,
label = '$ \\langle \\frac{d}{dx}V \\rangle - 2 \\gamma \\langle p \\rangle $' ,linewidth=1.5)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})
#ax.set_ylim(-0.8,0.8)
ax.set_xlabel('t')
#ax.set_ylabel(' ')
ax.grid();
In [16]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , np.gradient( instance.X2_average , instance.dt) , '-',
label='$\\frac{d}{dt} \\langle x^2 \\rangle$' , color = 'red', linewidth=1.5)
ax.plot( instance.timeRange , \
2*instance.XP_average/instance.mass, '--',label = '$\\frac{2}{m} \\langle xp \\rangle$',linewidth=1.5 )
#ax.set_xlim(0,3.5)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':24})
ax.set_xlabel('t')
ax.grid();
In [17]:
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot( instance.timeRange , np.gradient( instance.P2_average , instance.dt) , '-',
label = '$\\frac{d}{dt} \\langle p^2 \\rangle$',
color = 'red', linewidth=1.5)
ax.plot( instance.timeRange , \
-2*instance.PdPotentialdX_average\
+2.*instance.D_Theta \
-4*instance.gammaDamping*instance.P2_average
, '--',
label = '$- \\langle p\\frac{dV}{dx} +\\frac{dV}{dx} p \\rangle + 2 D_{\\theta}- 4\\gamma \\langle p^2 \\rangle $',
linewidth=1.5 )
#ax.set_xlim(-0.2,26)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(1.05, 1), loc=5, prop={'size':22})
ax.set_xlabel('t')
ax.grid();
In [18]:
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot( instance.timeRange , 2*np.gradient( instance.XP_average , instance.dt) ,
'-' ,label = '$\\frac{d}{dt} \\langle xp+px \\rangle$' , color = 'r' , linewidth=1.5 )
ax.plot( instance.timeRange , \
2*instance.P2_average/instance.mass \
-2*instance.XdPotentialdX_average \
-4*instance.gammaDamping*instance.XP_average
, '--' ,
label = '$\\frac{2}{m} \\langle p^2 \\rangle - 2 \\langle x \\frac{d}{dx}V \\rangle - 2 \\gamma \\langle xp + px \\rangle $'
,linewidth=1.5)
ax.legend(bbox_to_anchor=(1.05, 1), loc=5, prop={'size':22})
#ax.set_ylim(- 12 , 7)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [19]:
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot( instance.timeRange ,
np.sqrt(instance.X2_average - instance.X_average**2)*np.sqrt(instance.P2_average - instance.P_average**2)
, '-' , label = '$\\Delta x \\Delta p$' , linewidth=1.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=5, prop={'size':22})
#ax.set_ylim(-50, 0)
ax.set_xlabel('t')
#ax.set_xlim(0,41000)
ax.set_ylabel(' ')
ax.grid();
In [20]:
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot( instance.timeRange , instance.Hamiltonian_average
, '-' , label = '$Energy$' , linewidth=1.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=5, prop={'size':22})
#ax.set_ylim(3.48 , 3.52)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [21]:
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot( instance.timeRange , \
np.sqrt(instance.P2_average - instance.P_average**2) \
, '-' , label = '$\\langle p^2 \\rangle$',linewidth=2.)
#ax.plot( instance.timeRange , instance.X3_average - 2*gamma*instance.P_average , '-' ,
# label = '$-F-2\gamma <P>$' ,linewidth=2.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})
#ax.set_ylim(- 12 , 7)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [22]:
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot( instance.timeRange , \
np.sqrt(instance.X2_average - instance.X_average**2) \
, '-' , label = '$ \\langle x^2 \\rangle$',linewidth=2.)
#ax.plot( instance.timeRange , instance.X3_average - 2*gamma*instance.P_average , '-' ,
# label = '$-F-2\gamma <P>$' ,linewidth=2.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})
#ax.set_ylim(- 12 , 7)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [23]:
def dVdx(x):
return instance.dPotential(0,x)
In [24]:
xp_init = np.array([ instance.x_init, instance.p_init ])
trajectory = instance.SymplecticPropagator(
instance.dt, instance.timeSteps, xp_init, instance.gammaDamping ).T
In [25]:
xp_init = np.array([ instance.x_init, instance.p_init ])
trajectory_SmoothedP = instance.SymplecticPropagator_SmoothedP(
instance.dt, instance.timeSteps, xp_init, instance.gammaDamping ).T
In [26]:
x_min = -instance.X_amplitude
x_max = instance.X_amplitude - instance.dX
p_min = -instance.dP*instance.P_gridDIM/2.
p_max = instance.dP*instance.P_gridDIM/2. -instance.dP
print ' Quantum Ehrenfest trajectory vs classical trajectory'
print ' final time = ', instance.dt*instance.timeSteps
fig, ax = plt.subplots(figsize=(10, 10))
ax.plot( instance.X_average , instance.P_average , '-' ,color = 'g',
label = 'Fokker-Planck $(\\langle X \\rangle$, $\\langle P \\rangle )$ ', linewidth=1.5 )
#ax.plot( trajectory[0] , trajectory[1] , '--' , color='r',
# label = 'classical (X,P)', linewidth=1. )
ax.plot( trajectory_SmoothedP[0] , trajectory_SmoothedP[1] , '--' , color='b',
label = 'Newtonian $(X, P\, )$', linewidth=1. )
ax.set_xlabel('X')
ax.set_ylabel('P')
ax.set_aspect(1.5)
ax.set_xlim(-4.,4.)
ax.set_ylim(-4.,4.)
ax.contour(instance.Hamiltonian ,
np.arange(-45, 100, 2 ),origin='lower',extent=[x_min,x_max,p_min,p_max],
linewidths=0.25,colors='k')
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
#ax.grid();
ax.set_aspect(1.)
In [27]:
def NumpyPropertyOnly (obj):
"""
Class that will contain only picklable properties
"""
obj_ = dict()
for prop_name in dir(obj) :
prop = getattr(obj, prop_name)
if isinstance(prop, np.ndarray) : obj_[prop_name] = prop
return obj_
In [28]:
#pickle.dump( NumpyPropertyOnly(instance) , open( "X2.pickle", "wb" ) )
In [29]:
W = instance.WignerFunctionFromFile(1*instance.skipFrames)
In [30]:
instance.PlotWignerFrame( W.real ,
plotRange=((-6.,6) ,(-6,6)),
global_color=(-0.2, 0.2),
energy_Levels=(0, 20, 1), aspectRatio=1.);
In [31]:
quantumFileNameL3 = '/home/rcabrera/DATA/Rho2D/Ohmic/grid256_L3.0_epsilon0.5_amplitude20/'
quantumFileNameL3 += 'Rho2D_ODMDamping_XX_PInit3.0_Omega1_DTheta0.0143.hdf5'
In [32]:
quantumFileNameL6 = '/home/rcabrera/DATA/Rho2D/Ohmic/grid256_L6.0_epsilon0.5_amplitude20/'
quantumFileNameL6 += 'Rho2D_ODMDamping_XX_PInit3.0_Omega1_DTheta0.0143.hdf5'
In [33]:
fig, ax = plt.subplots(figsize=(10, 6))
quantum_timeRangeL3, quantum_X_averageL3 = instance.Ehrenfest_X_FromFile( fileName=quantumFileNameL3 )
quantum_timeRangeL6, quantum_X_averageL6 = instance.Ehrenfest_X_FromFile( fileName=quantumFileNameL6 )
ax.plot( quantum_timeRangeL3 , quantum_X_averageL3 , 'r-' ,
label='ODM L=3 $ \\langle x \\rangle$', linewidth=1.5 )
#ax.plot( quantum_timeRangeL6 , quantum_X_averageL6 , 'g-' ,
# label='ODM L=6 $\\langle x \\rangle$', linewidth=1.5 )
ax.plot( instance.timeRange , instance.X_average , '--' ,
label='Fokker-Planck $\\langle x \\rangle$', linewidth=1.5 )
axis_font = {'size':'24'}
ax.set_xlim(0,100)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(0.45, 1), loc=2, prop={'size':22})
ax.set_xlabel('$t \,\,(a.u.)$',**axis_font)
ax.set_ylabel('$\\langle x \\rangle \, (a.u.)$',**axis_font)
ax.grid();
fig.savefig('ODM-FokkerPlanckL3D0.143_X.pdf')
In [34]:
fig, ax = plt.subplots(figsize=(10, 6))
quantum_timeRangeL3, quantum_P_averageL3 = instance.Ehrenfest_P_FromFile( fileName=quantumFileNameL3 )
quantum_timeRangeL6, quantum_P_averageL6 = instance.Ehrenfest_P_FromFile( fileName=quantumFileNameL6 )
ax.plot( quantum_timeRangeL3 , quantum_P_averageL3 , 'r-' ,
label='ODM L=3 $\\langle p \\rangle$', linewidth=1.5 )
#ax.plot( quantum_timeRangeL6 , quantum_P_averageL6 , 'g-' ,
# label='ODM L=6 $\\langle p \\rangle$', linewidth=1.5 )
ax.plot( instance.timeRange , instance.P_average , '--' ,
label='Fokker-Planck $\\langle p \\rangle$', linewidth=1.5 )
axis_font = {'size':'24'}
ax.set_xlim(0,100)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(0.45, 1), loc=2, prop={'size':22})
ax.set_xlabel('$t \,\,(a.u.)$',**axis_font)
ax.set_ylabel('$\\langle p \\rangle \, (a.u.)$',**axis_font)
ax.grid();
fig.savefig('ODM-FokkerPlanckL3D0.143_P.pdf')
In [35]:
fig, ax = plt.subplots(figsize=(10, 6))
quantum_timeRangeL3, quantum_P2_averageL3 = instance.Ehrenfest_P2_FromFile( fileName=quantumFileNameL3 )
quantum_timeRangeL6, quantum_P2_averageL6 = instance.Ehrenfest_P2_FromFile( fileName=quantumFileNameL6 )
ax.plot( quantum_timeRangeL3 , quantum_P2_averageL3 , 'r-' ,
label='ODM L=3 $\\langle p^2 \\rangle$', linewidth=1.5 )
#ax.plot( quantum_timeRangeL6 , quantum_P2_averageL6 , 'g-' ,
# label='ODM L=6 $\\langle p^2 \\rangle$', linewidth=1.5 )
ax.plot( instance.timeRange , instance.P2_average , '--' ,
label='Fokker-Planck $\\langle p^2 \\rangle$', linewidth=1.5 )
axis_font = {'size':'24'}
ax.set_xlim(0,100)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(0.45, 1), loc=2, prop={'size':22})
ax.set_xlabel('$t \,\,(a.u.)$',**axis_font)
ax.set_ylabel('$\\langle p^2 \\rangle \, (a.u.)$',**axis_font)
ax.grid();
fig.savefig('ODM-FokkerPlanckL3D0.143_P2.pdf')
In [36]:
fig, ax = plt.subplots(figsize=(10, 6))
quantum_timeRangeL3, quantum_X2_averageL3 = instance.Ehrenfest_X2_FromFile( fileName=quantumFileNameL3 )
quantum_timeRangeL6, quantum_X2_averageL6 = instance.Ehrenfest_X2_FromFile( fileName=quantumFileNameL6 )
ax.plot( quantum_timeRangeL3 , quantum_X2_averageL3 , 'r-' ,
label='ODM L=3 $\\langle x^2 \\rangle$', linewidth=1.5 )
#ax.plot( quantum_timeRangeL6 , quantum_X2_averageL6 , 'g-' ,
# label='ODM L=6 $\\langle x^2 \\rangle$', linewidth=1.5 )
ax.plot( instance.timeRange , instance.X2_average , '--' ,
label='Fokker-Planck $\\langle x^2 \\rangle$', linewidth=1.5 )
axis_font = {'size':'24'}
ax.set_xlim(0,100)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(0.45, 1), loc=2, prop={'size':22})
ax.set_xlabel('$t \,\,(a.u.)$',**axis_font)
ax.set_ylabel('$\\langle x^2 \\rangle \, (a.u.)$',**axis_font)
ax.grid();
fig.savefig('ODM-FokkerPlanckL3D0.143_X2.pdf')
In [37]:
fig, ax = plt.subplots(figsize=(10, 6))
quantum_timeRangeL3, quantum_Hamiltonian_averageL3 = instance.Ehrenfest_Hamiltonian_FromFile( fileName=quantumFileNameL3 )
quantum_timeRangeL6, quantum_Hamiltonian_averageL6 = instance.Ehrenfest_Hamiltonian_FromFile( fileName=quantumFileNameL6 )
ax.plot( quantum_timeRangeL3 , quantum_Hamiltonian_averageL3 , 'r-' ,
label='ODM L=3 $\\langle H \\rangle$', linewidth=1.5 )
#ax.plot( quantum_timeRangeL6 , quantum_Hamiltonian_averageL6 , 'g-' ,
# label='ODM L=6 $\\langle H \\rangle$', linewidth=1.5 )
ax.plot( instance.timeRange , instance.Hamiltonian_average , '--' ,
label='Fokker-Planck $\\langle H \\rangle$', linewidth=1.5 )
ax.set_xlim(0,100)
#ax.set_ylim(-1.,1.2)
ax.legend(bbox_to_anchor=(0.45, 1), loc=2, prop={'size':22})
ax.set_xlabel('$t \,\,(a.u.)$',**axis_font)
ax.set_ylabel('$\\langle H \\rangle \, (a.u.)$',**axis_font)
ax.grid();
fig.savefig('ODM-FokkerPlanckL3D0.143_H.pdf')
In [38]:
W_end = instance.WignerFunctionFromFile( instance.timeSteps )
In [39]:
instance.PlotWignerFrame( W_end ,
plotRange=((-6.,6) ,(-6,6)),
global_color=(-0.2, 0.2),
energy_Levels=(0, 20, 1), aspectRatio=1.);
In [40]:
P_Op = instance.OperatorP_XBasis()
X_Op = instance.OperatorX_XBasis()
H_Op = P_Op.dot(P_Op)/(2*instance.mass) + 0.5*instance.mass*instance.omega**2*X_Op.dot(X_Op)
In [41]:
kT = instance.D_Theta/(2*instance.mass*instance.gammaDamping)
print ' kT = ', kT
In [42]:
BoltzmannRho = linalg.expm( - H_Op/kT )
norm = np.trace( BoltzmannRho )
BoltzmannRho /= norm
BoltzmannProbX = np.diagonal(BoltzmannRho).real
fig, ax = plt.subplots(figsize=(10, 5))
ProbX_end = np.sum( fftpack.fftshift(W_end.real) ,axis=0)*instance.dP
ax.semilogy( instance.X_range , ProbX_end ,label='Fokker-Planck: Probability' )
ax.semilogy( instance.X_range , BoltzmannProbX ,label='Quantum Boltzmann : Probability' )
ax.set_xlim( -5,5 )
ax.set_ylim( 1e-8, 2 )
ax.set_xlabel('x')
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[42]:
In [43]:
BoltzmanW = np.exp( -instance.Hamiltonian/kT )
norm = np.sum(BoltzmanW)*instance.dX*instance.dP
BoltzmanW /= norm
BoltzmannWProbX = np.sum( BoltzmanW ,axis=0 )*instance.dP
fig, ax = plt.subplots(figsize=(10, 5))
ProbX_end = np.sum( fftpack.fftshift(W_end.real) ,axis=0)*instance.dP
ax.semilogy( instance.X_range , ProbX_end ,label='Fokker-Planck: Probability' )
ax.semilogy( instance.X_range , BoltzmannWProbX ,label='Classical Boltzmann : Probability' )
ax.set_xlim( -5,5 )
ax.set_ylim( 1e-8, 2 )
ax.set_xlabel('x')
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[43]:
In [44]:
fig, ax = plt.subplots(figsize=(10, 5))
ProbP_end = np.sum( fftpack.fftshift(W_end.real) ,axis=1)*instance.dX
ax.semilogy( instance.X_range , ProbX_end ,label='Fokker-Planck: Probability ' )
ax.set_xlim( -5,5 )
ax.set_ylim( 1e-8, 2 )
ax.set_xlabel('p')
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[44]:
In [45]:
def WignerStationary( x,p , kT):
return np.exp( -( p**2/instance.mass + instance.mass*instance.omega**2*x**2 )/( instance.omega/np.tanh( instance.omega/2./kT ) ) )
In [46]:
W_stat = WignerStationary( instance.X , instance.P , kT)
In [47]:
instance.PlotWignerFrame( W_stat ,
plotRange=((-6.,6) ,(-6,6)),
global_color=(-0.2, 0.2),
energy_Levels=(0, 20, 1), aspectRatio=1.);
In [47]:
In [47]:
In [47]: