In [1]:
import pickle
import h5py
import numpy as np
import scipy.fftpack as fftpack
from scipy.special import hyp1f1
import pylab as plt
import scipy.linalg as linalg
from scipy import integrate
import time
#import ctypes
from string import Template
import sys
import os
import cv2
import pycuda.gpuarray as gpuarray
import pycuda.driver as cuda
from pycuda.compiler import SourceModule
import pycuda.autoinit
#import cublas
import scikits.cuda.linalg as cu_linalg
cu_linalg.init()
from pywignercuda_path import SetPyWignerCUDA_Path
SetPyWignerCUDA_Path()
from GPU_Rho2D_MUBasis import *
import matplotlib
%matplotlib inline
In [2]:
#Defining parameters
X_gridDIM = 256
X_amplitude = 20
dt= 0.005
timeSteps = 20000
skipFrames = 50
gammaRenormalization = 1. #15
print ' Evolution time ',timeSteps*dt
In [3]:
# Defining the instance
class Rho2D_Dissipation_X2(Rho2D_Dissipation):
"""
Instance for dissipative dynamics
"""
def __init__(self):
self.x_init = 0.
self.p_init = 0.
self.omega = 1
self.mass = 1
self.gammaDamping = 0.05
self.D_Theta = 0.0
self.L_material = 3.
self.epsilon = 0.5
epsilon = self.epsilon
#
constantX2 = 0.5*self.mass*self.omega**2
self.potentialString = '{0}*pow(x,2)'.format( constantX2 )
self.dPotentialString = '2.*{0}*x'.format( constantX2 )
dampingType = 'Ohmic'
# Data output file
self.fileName ='/home/rcabrera/DATA/Rho2D/'+dampingType
self.fileName +='/grid{2}_L{0}_epsilon{1}_amplitude{3}'.format( self.L_material,epsilon,X_gridDIM,X_amplitude)
self.fileName +='/Rho2D_ODMDamping_XX_PInit{0}_Omega{1}.hdf5'.format(self.p_init,self.omega)
# Lindbladian operator A
self.A_fileName = '/home/rcabrera/Documents/source/python/PyWignerCUDA/instances/Rho2D/'+dampingType
self.A_fileName += '/grid{2}_L{0}_epsilon{1}_amplitude{3}'.format(self.L_material,epsilon,X_gridDIM,X_amplitude)
self.A_fileName += '/Lindblad_A.hdf5'
Rho2D_Dissipation.__init__(
self,X_gridDIM,X_amplitude,dt,timeSteps,skipFrames)
# Setting initial state
n = 0
self.Rho_init = self.Rho_HarmonicOscillator(n,self.omega,self.x_init, self.p_init)
def f_Damping(self, P):
return P*P/np.sqrt( P*P + self.epsilon**2 )
def df_Damping(self, P):
return (P*P*P + 2*P*self.epsilon**2 )/( P*P + self.epsilon**2 )**(1.5)
def Ehrenfest_X2_QuantumCorrection(self,P):
hBar = self.hBar
L = self.L_material
gamma = self.gammaDamping
return 0.5*gamma*self.hBar*L*( P*P + 2*self.epsilon**2 )**2/( P*P + self.epsilon**2 )**(2.5)
In [4]:
instance = Rho2D_Dissipation_X2()
In [6]:
energySpectrum = np.sort( linalg.eigvals( instance.Hamiltonian ).real )
print ' dX = ',instance.dX , ' dP =',instance.dP
print ' X_Amplitude = ',instance.X_amplitude , ' P_Amplitude = ',instance.P_amplitude
print ' '
print ' Harm Osc Potential with energy spectrum '
fig, ax = plt.subplots(figsize=(10, 10))
X_range = instance.X_Range()
ax.plot( X_range , instance.Potential(0,X_range ) ,color = 'r')
map( lambda x: ax.plot( [-X_amplitude,X_amplitude], [x,x] , color='gray' ) ,energySpectrum)
ax.set_xlim(-10,10)
ax.set_ylim(0,30)
ax.set_xlabel('x')
ax.set_ylabel('V')
ax.set_aspect(1./2.)
#ax.grid('on')
In [7]:
fig, ax = plt.subplots(figsize=(10, 10))
ax.plot( instance.P_Range() , np.sign(instance.P_Range() )*instance.f_Damping( instance.P_Range() ) ,
'.', label=' sign(p) f(p) ' )
ax.plot( instance.P_Range() , instance.f_Damping( instance.P_Range() ) ,
'-', label=' f(p) ' )
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
ax.set_xlabel('$p$')
Out[7]:
In [8]:
instance.fileName
Out[8]:
In [9]:
if os.path.isfile( instance.A_fileName )==True:
print 'Loading A file'
A = instance.LoadLindbladianOperatorA( instance.A_fileName )
else:
print 'Starting the computation of Ohmic Lindbladian operator A'
initial_time = time.time()
A = instance.A_Damping( )
print "--- Calculation time ", time.time() - initial_time , ' s '
instance.SaveLindbladianOperatorA( instance.A_fileName , A )
In [10]:
initial_time = time.time()
instance.Run( A, gammaRenormalization )
print "--- Calculation time ", time.time() - initial_time , ' s '
In [11]:
print ' <X> = ', np.trace(instance.Rho_init.dot(instance.X)).real, ' -> ', np.trace(instance.Rho_end.dot(instance.X)).real
print ' <P> = ', np.round(np.trace(instance.Rho_init.dot(instance.P)).real,14) , ' -> ', np.trace(instance.Rho_end.dot(instance.P)).real
In [12]:
if os.path.isfile( instance.fileName )==True:
instance.Load_Ehrenfest()
In [13]:
print ' '
fig, ax = plt.subplots(figsize=(10, 4))
print ' normalization = ', np.trace( instance.Rho_end.real ), ' time = ', instance.timeRange[-1]
X_range = instance.X_Range()
maxProb = np.max( np.diagonal(instance.Rho_init).real/instance.dX )
ax.plot( X_range , np.log( np.diagonal(instance.Rho_init).real/instance.dX/maxProb ) , '--' ,label = 'init')
ax.plot( X_range , np.log( np.diagonal(instance.Rho_end).real/instance.dX/maxProb ) , '-' ,label = 'final')
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('Log Probability')
ax.set_ylim(-12,1)
#ax.set_xlim(-10,10)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[13]:
In [14]:
print ' '
fig, ax = plt.subplots(figsize=(10, 4))
print ' normalization = ', np.trace( instance.Rho_end.real ), ' time = ', instance.timeRange[-1]
P_range = instance.P_Range()
Rho_init_PRepresent = instance.X_To_P_Basis( instance.Rho_init )
Rho_end_PRepresent = instance.X_To_P_Basis( instance.Rho_end )
ax.semilogy( P_range , np.diagonal(Rho_init_PRepresent).real/instance.dP , '--' ,label = 'init')
ax.semilogy( P_range , np.diagonal(Rho_end_PRepresent).real/instance.dP , '-' ,label = 'final')
ax.grid()
ax.set_xlabel('p')
ax.set_ylabel('Probability Log')
#ax.set_xlim(-10,10)
ax.set_ylim(1e-6,1)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[14]:
In [15]:
print ' '
fig, ax = plt.subplots(figsize=(10, 4))
print ' normalization = ', np.trace( instance.Rho_end.real ), ' time = ', instance.timeRange[-1]
P_range = instance.P_Range()
Rho_init_PRepresent = instance.X_To_P_Basis( instance.Rho_init )
Rho_end_PRepresent = instance.X_To_P_Basis( instance.Rho_end )
ax.plot( P_range , np.diagonal(Rho_init_PRepresent).real/instance.dP , '--' ,label = 'init')
ax.plot( P_range , np.diagonal(Rho_end_PRepresent).real/instance.dP , '-' ,label = 'final')
ax.grid()
ax.set_xlabel('p')
ax.set_ylabel('Probability')
ax.set_xlim(-10,10)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[15]:
In [16]:
fig, ax = plt.subplots(figsize=(10, 4))
print ' normalization = ', np.trace( instance.Rho_end.real ),' time = ',instance.timeRange[-1]
Rho_init_XBasis = instance.Rho_init
Rho_end_XBasis = instance.Rho_end
ax.plot( instance.X_Range(), np.diagonal(Rho_init_XBasis).real/instance.dX ,label = 'init')
ax.plot( instance.X_Range(), np.diagonal(Rho_end_XBasis).real/instance.dX ,label = 'final')
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('Probability')
#ax.set_xlim(-15,20)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[16]:
In [17]:
def plot_Ehrenfest_dxdt():
time_index_max = 60000
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.timeRange , np.gradient( instance.X_average, instance.dt).real , '-',
label = '$\\frac{d}{dt} \\langle x \\rangle $' ,color = 'red', linewidth=2.5)
ax.plot( instance.timeRange , instance.P_average , '--' ,
label = '$\\frac{1}{m}\\langle p \\rangle$', linewidth=2.5 )
axis_font = {'fontname':'Times', 'size':'24'}
ax.set_xlabel(r'$t (a.u.)$',**axis_font)
#ax.set_xlim(0,3.5)
#ax.set_ylim(-1,1)
ax.legend(bbox_to_anchor=(0.7, 0.99), loc=2, prop={'size':22})
#ax.set_xlabel('t')
ax.grid();
return fig
In [18]:
plot_Ehrenfest_dxdt();
In [19]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , np.gradient( instance.P_average, instance.dt).real ,'-' ,
label = '$\\frac{d}{dt}p$' ,color = 'r' , linewidth=2.5)
ax.plot( instance.timeRange , -instance.dPotentialdX_average - 2*instance.gammaDamping*instance.P_average , '--' ,
label = '$-\\frac{d}{dx}V - 2 \\gamma sign(p) |f(p)| $' ,linewidth=2.)
axis_font = {'fontname':'Times', 'size':'24'}
ax.set_xlabel(r'$t (a.u.)$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':22})
ax.set_ylim(- 2.5 , 2.5)
#ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [20]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange, np.gradient( instance.X2_average, instance.dt).real ,'-' ,
label = '$\\frac{d}{dt}x^2$' ,color = 'r' , linewidth=2.)
ax.plot( instance.timeRange,
2*instance.XP_average/instance.mass + instance.delta_average,
'--' ,
label = '$\\frac{xp + px}{m} + \\frac{\\gamma\\hbar L}{2} \\frac{ \\left(\\frac{df(p)}{dp}\\right)^2 }{f(p)} $',
linewidth=2. )
axis_font = {'fontname':'Times', 'size':'24'}
ax.set_xlabel(r'$t (a.u.)$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':20} )
#ax.set_ylim(- 1 , 1)
#ax.set_xlim( 0 , 10)
#ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [21]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , np.gradient( instance.P2_average, instance.dt).real,
'-' , label = '$ \\frac{d}{dt} p^2 $',
color = 'r' , linewidth=2.)
ax.plot( instance.timeRange ,
-2.*instance.PdPotentialdX_average + \
instance.alpha_average
, '--' ,
label = '$-(p\\frac{dV}{dx} +\\frac{dV}{dx} p ) - 2\\gamma f( |p| )\\left( 2|p| - \\frac{\\hbar}{L} \\right) $',
linewidth=2.)
axis_font = {'fontname':'Times', 'size':'24'}
ax.set_xlabel(r'$t (a.u.)$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':20} )
#ax.set_ylim(- 0.1 , 0.1)
#ax.set_xlabel('t')
#ax.set_ylabel(' ')
ax.grid();
In [22]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , + 2*instance.gammaDamping*instance.f_Damping_average/instance.L_material
, '-' ,
label = '$ 2\\gamma \\hbar f(|p|)/L $',
linewidth=1.)
axis_font = {'fontname':'Times', 'size':'24'}
ax.set_xlabel(r'$t (a.u.)$',**axis_font)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, prop={'size':20} )
#ax.set_ylim(- 7 , 6)
#ax.set_xlabel('t')
#ax.set_ylabel(' ')
ax.grid();
In [23]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.timeRange , 2*np.gradient( instance.XP_average, instance.dt).real ,
'-' ,label = '$\\frac{d}{dt}(xp+px)$' ,
color = 'r' , linewidth=2.)
ax.plot( instance.timeRange , \
2*instance.P2_average/instance.mass \
-2*instance.XdPotentialdX_average
+instance.beta_average \
, '--' , label = '$\\frac{2}{m} p^2 - 2 x \\frac{d}{dx}V -4 \\gamma\, sign(p) f(|p|) x $',
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 [24]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.X2_average,
'-' ,label = '$<x^2>$' ,linewidth=2.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
#ax.set_ylim(- 6 , 5)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [25]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.P2_average,
'-' ,label = '$<p^2>$' ,linewidth=2.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
#ax.set_ylim(- 6 , 5)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [26]:
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot( instance.XP_average,
'-' ,label = '$<xp>$' ,linewidth=2.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
ax.set_ylim(- 1 , 1)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [27]:
fig, ax = plt.subplots(figsize=(10, 4))
sx = instance.X2_average - instance.X_average**2
sp = instance.P2_average - instance.P_average**2
ax.plot( np.sqrt(sx*sp) , '-' ,label = '$\\sigma_x^2 \\sigma^2_p$' ,linewidth=2.)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
#ax.set_ylim(0 , 2)
ax.set_xlabel('t')
ax.set_ylabel(' ')
ax.grid();
In [28]:
def EnergyWigner(n):
gridDIM = instance.gridDIM*2
Rho = instance.LoadRho(instance.fileName, n )
W = instance.DensityMatrix_To_Wigner(Rho)
X_amplitude = instance.X_amplitude*np.sqrt(2)
dX = instance.dX*np.sqrt(2)/2.
x_min = -X_amplitude
x_max = X_amplitude - dX
P_amplitude = instance.P_amplitude/np.sqrt(2)
dP = instance.dP/np.sqrt(2)/2.
p_min = -P_amplitude
p_max = P_amplitude - dP
X_Range = np.linspace( -X_amplitude, X_amplitude-dX, gridDIM )
P_Range = np.linspace( -P_amplitude, P_amplitude-dP, gridDIM )
X = X_Range[np.newaxis,:]
P = P_Range[:,np.newaxis]
H = P**2/(2*instance.mass) + instance.Potential(0,X)
return np.sum( W*H )*dX*dP
In [29]:
skipFactor = 10
energyW_list = []
for n in range( 0, instance.timeSteps, skipFactor*instance.skipFrames ):
print ' n = ',n
energyW_list.append( EnergyWigner(n) )
energyW_list = np.array( energyW_list )
energyW_timeLine = instance.dt*np.array( range( 0, instance.timeSteps, skipFactor*instance.skipFrames ) )
In [30]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(instance.timeRange, instance.P2_average/(2*instance.mass) + instance.Potential_average,
color='0.5', label = 'Total Energy (Rho)')
ax.plot(energyW_timeLine, energyW_list.real, '.',
color='red', label = 'Total Energy (Wigner)')
ax.set_ylim(0 , 5)
ax.set_xlabel('t')
ax.grid();
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[30]:
In [31]:
def Purity(n):
Rho = instance.LoadRho(instance.fileName, n )
return np.trace( Rho.dot(Rho) ).real
In [32]:
purity_list = []
for n in range( 0, instance.timeSteps, instance.skipFrames ):
purity_list.append( Purity(n) )
purity_list = np.array(purity_list )
purity_timeLine = instance.dt*np.array( range( 0, instance.timeSteps, instance.skipFrames ) )
In [33]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( purity_timeLine, purity_list ,label='Purity' )
ax.set_xlabel('t')
ax.set_ylabel('Purity')
ax.grid();
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[33]:
In [34]:
def dVdx(x):
return instance.dPotential(0,x)
In [35]:
# Approximate solution disregarding epsilon
trajectory = instance.SymplecticPropagator(
instance.mass, instance.gammaDamping, instance.dt, instance.timeSteps, dVdx,
np.array([ instance.x_init, instance.p_init ]) ).T
In [36]:
# python ODE solution
def dxdt_dpdt( Y , t):
return [ Y[1]/instance.mass , -dVdx(Y[0]) - instance.gammaDamping*Y[1]*np.abs(Y[1])/np.sqrt( Y[1]**2 + instance.epsilon ) ]
timeRange = np.arange(0, instance.dt*instance.timeSteps , instance.dt)
trajectory2 = integrate.odeint( dxdt_dpdt , [instance.x_init, instance.p_init], timeRange ).T
In [37]:
print 'final time = ' , instance.timeRange[-1]
p_min = -instance.dP *instance.gridDIM/2.
p_max = instance.dP *instance.gridDIM/2.
x_min = -instance.X_amplitude
x_max = instance.X_amplitude
x_range = instance.X_Range()
p_range = instance.P_Range()
X = x_range[np.newaxis,:]
P = p_range[:,np.newaxis]
H = P**2/(2*instance.mass) + instance.Potential(0,X)
fig, ax = plt.subplots(figsize=(10, 10))
ax.plot( instance.X_average , instance.P_average , '--' , label = 'quantum(<X>,<P>)', linewidth=2 )
#ax.plot( trajectory[0] , trajectory[1] , '.' , color='g', label = 'classical (X,P)', linewidth=1 )
ax.plot( trajectory[0] , trajectory[1] , '-' , color='r',
label = 'classical (X,P)', linewidth=1 )
ax.contour(H, np.arange(0, 4, 0.5 ),origin='lower',extent=[x_min,x_max,p_min,p_max],
linewidths=0.25,colors='k')
ax.set_xlabel('X')
ax.set_ylabel('P')
ax.set_aspect(1.)
ax.set_xlim(-3.,3.)
ax.set_ylim(-3.,3.)
ax.grid();
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[37]:
In [38]:
def PlotWigner( W2 , messages=True,saveFigure=False):
global_min = -0.08
global_max = 0.08
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, 512)
#wigner_cmap = LinearSegmentedColormap('wigner_colormap', wigner_cdict, 512)
plt.clf()
fig, ax = plt.subplots(figsize=(10, 8))
dp = instance.dP
dx = instance.dX
p_min = -instance.P_amplitude + dp/2.
p_max = instance.P_amplitude - dp/2.
x_min = -instance.X_amplitude + dx/2.
x_max = instance.X_amplitude - dx/2.
norm = np.sum(W2)*dp*dx
W2 /= norm
cax = ax.imshow( W2.real ,
extent = [ x_min , x_max , p_min , p_max ],
vmin= global_min, vmax=global_max,
origin='lower',interpolation='none',cmap=wigner_cmap)
X = instance.X_Range()[np.newaxis,:]
P = instance.P_Range()[:,np.newaxis]
H = P**2/(2*instance.mass) + instance.Potential(0,X)
ax.contour( H , energySpectrum , origin='lower',extent=[x_min,x_max,p_min,p_max],
linewidths=0.25,colors='k')
ax.set_xlim( -5,5 )
ax.set_ylim( -5,5 )
ax.set_xlabel('x')
ax.set_ylabel('p')
ax.grid('on')
ax.set_aspect(1)
cbar = fig.colorbar( cax , ticks=[-0.3, -0.2,-0.1, 0, 0.1, 0.2 , 0.3])
if saveFigure != False:
fig.savefig( saveFigure )
if messages == True:
print ' max W = ', np.max(W2.real)
print ' min W = ', np.min(W2.real)
print ' norm = ', np.sum(W2).real*dp*dx
print ' Energy =', np.sum( H*W2 ).real*dx*dp
#return fig
In [39]:
def PlotWignerScaled( W2 , messages=True,saveFigure=False):
"""
This functions must be called when plotting Wigner functions are calculated from the density matrix
"""
global_min = -0.08
global_max = 0.08
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, 512)
plt.clf()
fig, ax = plt.subplots( figsize=(10, 8) )
gridDIM = W2.shape[0]
X_amplitude = instance.X_amplitude*np.sqrt(2)
dX = instance.dX*np.sqrt(2)/2.
x_min = -X_amplitude
x_max = X_amplitude - dX
P_amplitude = instance.P_amplitude/np.sqrt(2)
dP = instance.dP/np.sqrt(2)/2.
p_min = -P_amplitude
p_max = P_amplitude - dP
X_Range = np.linspace( -X_amplitude, X_amplitude-dX, gridDIM )
P_Range = np.linspace( -P_amplitude, P_amplitude-dP, gridDIM )
X = X_Range[np.newaxis,:]
P = P_Range[:,np.newaxis]
cax = ax.imshow( W2.real ,
extent = [ x_min , x_max , p_min , p_max ],
vmin= global_min, vmax=global_max,
origin='lower',interpolation='none',cmap=wigner_cmap)
H = P**2/(2*instance.mass) + instance.Potential(0,X)
CS = ax.contour( H , energySpectrum[0:12] , origin='lower',extent=[x_min,x_max,p_min,p_max],
linewidths=0.5,colors='green')
ax.clabel( CS, inline=1,fontsize=12, fmt='%1.2f')
ax.set_xlim( -5,5 )
ax.set_ylim( -5,5 )
ax.set_xlabel('x')
ax.set_ylabel('p')
ax.grid('on')
ax.set_aspect(1)
cbar = fig.colorbar( cax , ticks=[-0.3, -0.2,-0.1, 0, 0.1, 0.2 , 0.3])
if saveFigure != False:
fig.savefig( saveFigure )
if messages == True:
print ' max W = ', np.max(W2.real)
print ' min W = ', np.min(W2.real)
print ' norm = ', np.sum(W2).real*dP*dX
print ' Energy (Wigner) =', np.sum( H*W2 ).real*dX*dP
print ' <xp> (Wigner) =', np.sum( W2*X*P).real*dX*dP
In [40]:
W0 = instance.DensityMatrix_To_Wigner( instance.LoadRho(instance.fileName,0) )
In [41]:
fig0 = PlotWignerScaled(W0);
In [42]:
W10000 = instance.DensityMatrix_To_Wigner( instance.LoadRho(instance.fileName, instance.timeSteps ) )
In [43]:
PlotWignerScaled(W10000);
In [44]:
WW = instance.DensityMatrix_To_Wigner(
instance.Rho_HarmonicOscillator( 2, instance.omega, 0, 0)
).real
In [45]:
PlotWignerScaled(WW)
In [46]:
PlotWigner(
instance.Wigner_HarmonicOscillator( 2 ,instance.omega,0,0)
);
In [47]:
nMax=60
projectionsOdd = []
for n in range(1,nMax,2):
RhoHO = instance.Rho_HarmonicOscillator(n , instance.omega , 0 , 0 )
Rho = instance.LoadRho(instance.fileName, instance.timeSteps)
projectionsOdd.append( np.trace( Rho.dot(RhoHO) ) )
projectionsOdd = np.array(projectionsOdd)
In [48]:
projectionsEven = []
for n in range(0,nMax,2):
RhoHO = instance.Rho_HarmonicOscillator(n , instance.omega, 0 , 0 )
Rho = instance.LoadRho(instance.fileName, instance.timeSteps)
projectionsEven.append( np.trace( Rho.dot(RhoHO) ) )
projectionsEven = np.array(projectionsEven)
In [49]:
plt.semilogy( range(1,nMax,2), projectionsOdd.real ,'.', label='Odd states')
plt.semilogy( range(0,nMax,2), projectionsEven.real ,'r.',label='Even states')
plt.legend(bbox_to_anchor=(1.05, 1), loc=2)
plt.xlabel('n')
Out[49]:
In [50]:
RhoSum = np.zeros( (instance.gridDIM,instance.gridDIM) , dtype=np.complex128)
for n in range(0,nMax/2):
RhoSum += projectionsEven[n]*instance.Rho_HarmonicOscillator(2*n , instance.omega , 0 , 0 )
for n in range(1,nMax/2):
RhoSum += projectionsOdd[n]*instance.Rho_HarmonicOscillator(2*n+1 , instance.omega , 0 , 0 )
In [51]:
WSum = instance.DensityMatrix_To_Wigner( RhoSum )
In [52]:
PlotWignerScaled(WSum);
In [53]:
# avconv -r 10 -i fig%05d.png -b:v 1000k test.mp4
"""
m=0
for n in range(0,instance.timeSteps,instance.skipFrames):
print ' n = ',n
m+=1
W = instance.DensityMatrix_To_Wigner( instance.LoadRho(instance.fileName,n) )
PlotWigner( W , messages=False, saveFigure='/home/rcabrera/DATA/Rho2D/Ohmic/movie1/fig%05d'%m );
""";
In [53]:
In [ ]: