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 matplotlib.cm as cm
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 = 40000
skipFrames = 100
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 = 3.
self.omega = 1
self.mass = 1
self.gammaDamping = 0.07
self.D_Theta = 0.02
#self.kT = 0.1
#self.D_Theta = self.D_Function( self.kT )
self.L_material = 2.
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}_DTheta{2}.hdf5'.format(
self.p_init, self.omega, self.D_Theta )
# 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 D_Function(self,kT,Omega=None):
if Omega is None:
Omega = 0.1*self.omega
if Omega == 0:
return 2*self.gammaDamping*self.mass*kT
else:
return self.gammaDamping*self.mass*Omega/2./np.tanh( Omega/(4* kT) )
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 [5]:
fig, ax = plt.subplots(figsize=(8, 8))
kT_Range = np.linspace( 1e-6, 0.2 ,100 )
ax.plot( kT_Range , instance.D_Function( kT_Range, Omega=0.1 ) ,label='$\\Omega=0.1$')
ax.plot( kT_Range , instance.D_Function( kT_Range, Omega=0.2 ) ,label='$\\Omega=0.2$')
ax.plot( kT_Range , instance.D_Function( kT_Range, Omega=0.3 ) ,label='$\\Omega=0.3$')
ax.plot( kT_Range , instance.D_Function( kT_Range, Omega=0.4 ) ,label='$\\Omega=0.4$')
ax.set_xlabel('kT')
ax.set_ylabel('D')
ax.grid('on')
ax.legend(bbox_to_anchor=(1.05, 1), loc=2)
Out[5]:
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=1.5)
ax.plot( instance.timeRange , instance.P_average , '--' ,
label = '$\\frac{1}{m}\\langle p \\rangle$', linewidth=1.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=1.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=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':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=1.)
ax.plot( instance.timeRange,
2*instance.XP_average/instance.mass + instance.delta_average + 0*instance.D_Theta,
'--' ,
label = '$\\frac{xp + px}{m} + \\frac{\\gamma\\hbar L}{2} \\frac{ \\left(\\frac{df(p)}{dp}\\right)^2 }{f(p)} $',
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(- 1 , 1)
ax.set_xlim( 0 , 100)
#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=1.)
ax.plot( instance.timeRange ,
-2.*instance.PdPotentialdX_average + \
instance.alpha_average + 2*instance.D_Theta
, '--' ,
label = '$-(p\\frac{dV}{dx} +\\frac{dV}{dx} p ) - 2\\gamma f( |p| )\\left( 2|p| - \\frac{\\hbar}{L} \\right) + 2D$',
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(- 1 , 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=1.)
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=1.)
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 + dX/2.
x_max = X_amplitude - dX/2.
P_amplitude = instance.P_amplitude/np.sqrt(2)
dP = instance.dP/np.sqrt(2)/2.
p_min = -P_amplitude + dP/2.
p_max = P_amplitude - dP/2.
X_Range = np.linspace( -X_amplitude + dX/2., X_amplitude - dX/2., gridDIM )
P_Range = np.linspace( -P_amplitude + dP/2., P_amplitude - dP/2., 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]:
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 + dX/2.
x_max = X_amplitude - dX/2.
P_amplitude = instance.P_amplitude/np.sqrt(2)
dP = instance.dP/np.sqrt(2)/2.
p_min = -P_amplitude + dP/2.
p_max = P_amplitude - dP/2.
X_Range = np.linspace( -X_amplitude + dX/2., X_amplitude - dX/2., gridDIM )
P_Range = np.linspace( -P_amplitude + dP/2., P_amplitude - dP/2., 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 [41]:
def PlotWignerScaled_Hot( W2 , messages=True,saveFigure=False):
"""
This functions must be called when plotting Wigner functions are calculated from the density matrix
"""
global_min = -1e-6
global_max = 0.13
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 + dX/2.
x_max = X_amplitude - dX/2.
P_amplitude = instance.P_amplitude/np.sqrt(2)
dP = instance.dP/np.sqrt(2)/2.
p_min = -P_amplitude + dP/2.
p_max = P_amplitude - dP/2.
X_Range = np.linspace( -X_amplitude + dX/2., X_amplitude - dX/2., gridDIM )
P_Range = np.linspace( -P_amplitude + dP/2., P_amplitude - dP/2., 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)
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=cm.hot_r)
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')
axis_font = {'size':'34'}
#ax.clabel( CS, inline=1,fontsize=12, fmt='%1.2f')
ax.set_xlim( -5,5 )
ax.set_ylim( -5,5 )
ax.set_xlabel('x (a.u.)',**axis_font)
ax.set_ylabel('p (a.u.)',**axis_font)
plt.xticks(fontsize = 28)
plt.yticks(fontsize = 28)
#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
return fig
In [42]:
W0 = instance.DensityMatrix_To_Wigner( instance.LoadRho(instance.fileName,0) )
In [43]:
fig0 = PlotWignerScaled_Hot(W0);
In [44]:
W10000 = instance.DensityMatrix_To_Wigner( instance.LoadRho(instance.fileName, instance.timeSteps ) )
In [45]:
fig1 = PlotWignerScaled_Hot(W10000);
In [46]:
instance.D_Theta
Out[46]:
In [47]:
fig1.savefig('finalState_L3_D0.01429_P3.pdf')
$D = 2 m k T \gamma$
In [48]:
effective_kT = 1.2908 #1.166
print '************* Effective kT = ',effective_kT,' *******************'
BoltzmannRho = linalg.expm( -instance.Hamiltonian/effective_kT )
norm = np.trace( BoltzmannRho )
BoltzmannRho /= norm
In [49]:
W_Boltzmann = instance.DensityMatrix_To_Wigner( BoltzmannRho )
In [50]:
fig2 = PlotWignerScaled_Hot(W_Boltzmann);
In [51]:
#fig2.savefig('BoltzmannState_kT1.166.pdf')
In [52]:
BoltzmannProb_x = np.diagonal(BoltzmannRho).real
print ' Boltzmann Norm = ', np.trace( BoltzmannRho ).real
print ' Purity Boltzmann = ', np.trace( BoltzmannRho.dot(BoltzmannRho) ).real
Rho_End = instance.LoadRho(instance.fileName, instance.timeSteps )
RhoEnd_Prob_x = np.diagonal( Rho_End ).real
print ' '
print ' Rho Norm = ', np.trace( Rho_End ).real
print ' Purity Rho_End = ', np.trace( Rho_End.dot(Rho_End) ).real
print ' '
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.X_Range() , BoltzmannProb_x ,label='Boltzmann Prob x' )
ax.plot( instance.X_Range() , RhoEnd_Prob_x ,label='Rho End Prob x' )
ax.set_xlim( -5,5 )
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('Prob')
Out[52]:
In [53]:
BoltzmannProb_p = np.diagonal( instance.X_To_P_Basis( BoltzmannRho ) ).real
print ' Boltzmann Norm = ', np.trace( BoltzmannRho ).real
Rho_End = instance.X_To_P_Basis( instance.LoadRho(instance.fileName, instance.timeSteps ) )
RhoEnd_Prob_p = np.diagonal( Rho_End ).real
print ' Rho Norm = ', np.trace( Rho_End ).real
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.X_Range() , BoltzmannProb_p ,label='Boltzmann Prob p' )
ax.plot( instance.X_Range() , RhoEnd_Prob_p ,label='Rho End Prob p' )
ax.set_xlim( -5,5 )
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
ax.grid()
ax.set_xlabel('p')
ax.set_ylabel('Prob')
Out[53]:
In [54]:
BoltzmannProb_x = np.diagonal(BoltzmannRho).real
print ' Boltzmann Norm = ', np.trace( BoltzmannRho ).real
print ' Purity Boltzmann = ', np.trace( BoltzmannRho.dot(BoltzmannRho) ).real
Rho_End = instance.LoadRho(instance.fileName, instance.timeSteps )
RhoEnd_Prob_x = np.diagonal( Rho_End ).real
print ' '
print ' Rho Norm = ', np.trace( Rho_End ).real
print ' Purity Rho_End = ', np.trace( Rho_End.dot(Rho_End) ).real
print ' '
fig, ax = plt.subplots(figsize=(10, 5))
ax.semilogy( instance.X_Range() , BoltzmannProb_x ,label='Boltzmann Prob x' )
ax.semilogy( instance.X_Range() , RhoEnd_Prob_x ,label='Rho End Prob x' )
ax.set_xlim( -15,15 )
ax.set_ylim( 1e-6,2 )
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('Prob')
Out[54]:
In [55]:
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 [56]:
from scipy.optimize import curve_fit
def FindkT(Prob_p, P_Range, omega, mass):
"""
Function to find kT by maximazing the overlap with the Boltzmann state
"""
def BoltzmannPMarginal(P, kT):
if kT == 0:
C = 1/(omega*mass)
else:
C = np.tanh(0.5*omega/kT) / (omega*mass)
B = np.exp(-C*P**2)
B /= B.sum()
return B
return curve_fit(BoltzmannPMarginal, P_Range, Prob_p)[0][0]
# test
FindkT(np.diagonal( Rho_End ).real, instance.P_Range(), instance.omega , instance.mass)
Out[56]:
In [57]:
#dBoltzmannRho = BoltzmannRho
#for i in range(50):
# dBoltzmannRho = linalg.sqrtm( dBoltzmannRho )
In [58]:
# Testing Boltzmann decomposition
# np.abs( BoltzmannRho - BoltzmannRho_Cho.dot(BoltzmannRho_Cho.conj().T) ).max()
#V = dBoltzmannRho.dot( linalg.expm( -1j*instance.Hamiltonian*instance.dt ) )
In [59]:
#RhoInit2 = instance.LoadRho(instance.fileName, 10*instance.skipFrames )
#Rho_Cho = linalg.cholesky( RhoInit2 ).conj().T
In [60]:
# Testing Initial state decomposition
#np.abs( RhoInit2 - B_init.dot( B_init.conj().T ) ).max()
In [61]:
"""x_average = []
p_average = []
for t_index in range(1000):
norm = np.trace( RhoInit2 )
RhoInit2 /= norm
x_average.append( np.sum( RhoInit2*instance.X ) )
p_average.append( np.sum( RhoInit2*instance.P ) )
RhoInit2 = V.dot( RhoInit2.dot( V.conj().T ) )
x_average = np.array(x_average).real
p_average = np.array(p_average).real """;
In [72]:
RhoInit = instance.Rho_init + (1e-2)*BoltzmannRho
In [98]:
instance.Run_To_Boltzmann( RhoInit , effective_kT )
Out[98]:
In [99]:
WTerminal = instance.DensityMatrix_To_Wigner( instance.RhoTerminal_end )
In [100]:
PlotWignerScaled_Hot( WTerminal );
In [101]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.X_Range() , BoltzmannProb_p ,label='Boltzmann Prob p' )
ax.plot( instance.X_Range() , np.diagonal( instance.RhoTerminal_end ).real ,label='Terminal' )
ax.set_xlim( -5,5 )
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[101]:
In [102]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.X_average_ , 'r-' ,
label = ' Terminal $\\langle x \\rangle$', linewidth=1.5 )
ax.plot( instance.X_average , '--' ,
label = 'ODM $\\langle x \\rangle$', linewidth=1.5 )
#ax.set_xlim(0, 10000)
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[102]:
In [103]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.P_average_ , 'r-' ,
label = ' Terminal $\\langle x \\rangle$', linewidth=1.5 )
ax.plot( instance.P_average , '--' ,
label = 'ODM $\\langle x \\rangle$', linewidth=1.5 )
#ax.set_xlim(0, 100)
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[103]:
In [104]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.X2_average_ , 'r-' ,
label = ' Terminal $\\langle x \\rangle$', linewidth=1.5 )
ax.plot( instance.X2_average , '--' ,
label = 'ODM $\\langle x \\rangle$', linewidth=1.5 )
#ax.set_xlim(0, 100)
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[104]:
In [105]:
fig, ax = plt.subplots(figsize=(10, 5))
ax.plot( instance.P2_average_ , 'r-' ,
label = ' Terminal $\\langle x \\rangle$', linewidth=1.5 )
ax.plot( instance.P2_average , '--' ,
label = 'ODM $\\langle x \\rangle$', linewidth=1.5 )
#ax.set_xlim(0, 100)
ax.legend(bbox_to_anchor=(0.8, 1), loc=2)
Out[105]:
In [66]:
WW = instance.DensityMatrix_To_Wigner(
instance.Rho_HarmonicOscillator( 2, instance.omega, 0, 0)
).real
In [67]:
PlotWignerScaled(WW)
In [68]:
PlotWigner(
instance.Wigner_HarmonicOscillator( 2 ,instance.omega,0,0)
);
In [69]:
# 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 );
""";