1
#This file is part of QuTIP.
3
# QuTIP is free software: you can redistribute it and/or modify
4
# it under the terms of the GNU General Public License as published by
5
# the Free Software Foundation, either version 3 of the License, or
6
# (at your option) any later version.
8
# QuTIP is distributed in the hope that it will be useful,
9
# but WITHOUT ANY WARRANTY; without even the implied warranty of
10
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11
# GNU General Public License for more details.
13
# You should have received a copy of the GNU General Public License
14
# along with QuTIP. If not, see <http://www.gnu.org/licenses/>.
16
# Copyright (C) 2011, Paul D. Nation & Robert J. Johansson
18
###########################################################################
21
from scipy.integrate import *
25
# ------------------------------------------------------------------------------
26
# Wave function evolution using a ODE solver (unitary quantum evolution)
28
def wf_ode_evolve(tlist, H, psi0, expt_op_list):
30
@brief Evolve the wave function using an ODE solver
33
n_expt_op = len(expt_op_list)
34
expt_list = zeros([n_expt_op, n_tsteps])
35
dt = tlist[1]-tlist[0]
37
r = scipy.integrate.ode(mcwf_ode_func)
39
initial_vector = psi0.full()
40
#print "initial vector = ", initial_vector
42
r.set_integrator('zvode').set_initial_value(initial_vector, 0.0).set_f_params(H)
45
for m in range(0, n_expt_op):
46
expt_list[m,t_idx] += expect(expt_op_list[m], psi0)
48
while r.successful() and r.t < tlist[-1]:
54
# calculate all the expectation values: contribution from single trajectory
55
for m in range(0, n_expt_op):
56
expt_list[m,t_idx] += expect(expt_op_list[m], psi)
65
def mcwf_ode_func(t, psi, H):
66
dpsi = -1j * (H.data * psi)
70
# ------------------------------------------------------------------------------
71
# Master equation solver
73
def me_ode_evolve(tlist, H, rho0, c_op_list, expt_op_list):
75
@brief Evolve the denstiy matrix using an ODE solver
78
n_expt_op = len(expt_op_list)
80
expt_list = zeros([n_expt_op, n_tsteps], dtype=complex)
81
dt = tlist[1]-tlist[0]
84
# Got a wave function as initial state: convert to density matrix.
85
rho0 = rho0 * trans(rho0)
88
# construct liouvillian
90
L = -1j*(spre(H) - spost(H))
91
for m in range(0, n_op):
92
cdc = c_op_list[m].dag() * c_op_list[m]
93
L += spre(c_op_list[m])*spost(c_op_list[m].dag())-0.5*spre(cdc)-0.5*spost(cdc)
98
initial_vector = rho0.full().reshape(prod(rho0.shape),1)
99
r = scipy.integrate.ode(me_ode_func).set_integrator('zvode').set_initial_value(initial_vector, tlist[0]).set_f_params(L.data)
101
#r.set_integrator('zvode').set_initial_value(initial_vector, 0.0).set_f_params(L)
109
# while t_idx <= n_tsteps:
111
if not r.successful():
114
rho.data = r.y.reshape(rho0.shape)
116
# calculate all the expectation values
117
for m in range(0, n_expt_op):
118
expt_list[m,t_idx] = expect(expt_op_list[m], rho)
120
r.integrate(r.t + dt)
127
# evaluate drho(t)/dt according to the master eqaution
129
def me_ode_func(t, rho, L):
added
removed
qutip/wf_ode_evolve.py
