From 7dbc578435ba843eeccb862a2e3862153c42a818 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Man=C4=8Dal?= Date: Thu, 11 Jan 2024 15:41:18 +0100 Subject: [PATCH] First steps towards implementing state vector open quantum systems description --- quantarhei/__init__.py | 2 + quantarhei/qm/__init__.py | 2 + quantarhei/qm/hilbertspace/oqsstatevector.py | 38 ++++- .../liouvillespace/rates/tdredfieldrates.py | 2 +- quantarhei/qm/propagators/oqssvevolution.py | 55 ++++++++ quantarhei/qm/propagators/oqssvpropagator.py | 132 ++++++++++++++++++ .../qm/hilbertspace/oqsstatevector_test.py | 22 +++ 7 files changed, 246 insertions(+), 7 deletions(-) create mode 100644 quantarhei/qm/propagators/oqssvevolution.py create mode 100644 quantarhei/qm/propagators/oqssvpropagator.py diff --git a/quantarhei/__init__.py b/quantarhei/__init__.py index f81c484..7c44fb3 100644 --- a/quantarhei/__init__.py +++ b/quantarhei/__init__.py @@ -326,12 +326,14 @@ # from .qm.propagators.poppropagator import PopulationPropagator from .qm.propagators.svpropagator import StateVectorPropagator +from .qm import OQSStateVectorPropagator from .qm import ReducedDensityMatrixPropagator # # Evolutions (time-dependent operators) # from .qm.propagators.statevectorevolution import StateVectorEvolution +from .qm import OQSStateVectorEvolution from .qm import DensityMatrixEvolution from .qm import ReducedDensityMatrixEvolution diff --git a/quantarhei/qm/__init__.py b/quantarhei/qm/__init__.py index 7959bb9..5192d5c 100644 --- a/quantarhei/qm/__init__.py +++ b/quantarhei/qm/__init__.py @@ -80,11 +80,13 @@ # from .propagators.rdmpropagator import ReducedDensityMatrixPropagator from .propagators.svpropagator import StateVectorPropagator +from .propagators.oqssvpropagator import OQSStateVectorPropagator # # TIME EVOLUTIONS # from .propagators.statevectorevolution import StateVectorEvolution +from .propagators.oqssvevolution import OQSStateVectorEvolution from .propagators.dmevolution import DensityMatrixEvolution from .propagators.dmevolution import ReducedDensityMatrixEvolution diff --git a/quantarhei/qm/hilbertspace/oqsstatevector.py b/quantarhei/qm/hilbertspace/oqsstatevector.py index 552957a..184d14f 100644 --- a/quantarhei/qm/hilbertspace/oqsstatevector.py +++ b/quantarhei/qm/hilbertspace/oqsstatevector.py @@ -8,6 +8,7 @@ import numpy from ... import REAL, COMPLEX +from .operators import ReducedDensityMatrix class OQSStateVector(): """Represents a quantum mechanical state of an open system @@ -40,7 +41,7 @@ class OQSStateVector(): def __init__(self, dim=None, data=None): - self._initialited = False + self._initialized = False if data is not None: @@ -49,18 +50,19 @@ def __init__(self, dim=None, data=None): if len(ddat.shape) > 1: raise Exception("Data has to be a vector") - if dim != ddat.shape[0]: - print("Dimension specification differes from data: ignored.") + if dim is not None: + if dim != ddat.shape[0]: + print("Dimension specification different from data: ignored.") self.data = ddat self.dim = ddat.shape[0] - self._initialited = True + self._initialized = True elif dim is not None: self.dim = dim self.data = numpy.zeros(self.dim, dtype=REAL) - self._initialited = True + self._initialized = True def norm(self): @@ -68,7 +70,7 @@ def norm(self): """ - return numpy.dot(self.data,self.data) + return numpy.sqrt(numpy.dot(self.data,self.data)) def puredot(self, psi): @@ -77,5 +79,29 @@ def puredot(self, psi): """ return numpy.dot(self.data, psi.data) + + + def get_ReducedDensityMatrix(self, decoherence=False): + """ Converts the state vector into the density matrix + + Parameters + ---------- + decoherence : bool, optional + Should the density matrix contain decoherence?. The default is False. + + Returns + ------- + ReducedDensityMatrix object with real values. + Effectively it is in interaction picture. + + """ + data = numpy.zeros((self.dim, self.dim), dtype=REAL) + rho = ReducedDensityMatrix(data=data) + for ii in range(self.dim): + for jj in range(self.dim): + rho.data[ii,jj] = self.data[ii]*self.data[jj] + + return rho + \ No newline at end of file diff --git a/quantarhei/qm/liouvillespace/rates/tdredfieldrates.py b/quantarhei/qm/liouvillespace/rates/tdredfieldrates.py index 5300ea3..6811256 100644 --- a/quantarhei/qm/liouvillespace/rates/tdredfieldrates.py +++ b/quantarhei/qm/liouvillespace/rates/tdredfieldrates.py @@ -55,7 +55,7 @@ def __init__(self, ham, sbi, initialize=True, cutoff_time=None): raise Exception("First argument must be a Hamiltonian") if not isinstance(sbi,SystemBathInteraction): - raise Exception + raise Exception("Second argument must be a SystemBathInteraction") self._is_initialized = False self._has_cutoff_time = False diff --git a/quantarhei/qm/propagators/oqssvevolution.py b/quantarhei/qm/propagators/oqssvevolution.py new file mode 100644 index 0000000..bad4ef8 --- /dev/null +++ b/quantarhei/qm/propagators/oqssvevolution.py @@ -0,0 +1,55 @@ +# -*- coding: utf-8 -*- +import numpy + +from ... import REAL, COMPLEX +from ..propagators.dmevolution import ReducedDensityMatrixEvolution + +class OQSStateVectorEvolution: + + + def __init__(self, timeaxis=None, psii=None): + + + if timeaxis is not None: + + self.TimeAxis = timeaxis + + if psii is not None: + self.dim = psii.data.shape[0] + + self.set_initial_condition(psii) + else: + self.dim = 0 + + + + def set_initial_condition(self, psii): + """ + + + """ + self.dim = psii.data.shape[0] + self.data = numpy.zeros((self.TimeAxis.length, \ + self.dim),dtype=numpy.float64) + self.data[0,:] = psii.data + + + def get_ReducedDensityMatrixEvolution(self, decoherence=False): + """ Returns the corresponding reduced density matrix + + + The present implementation cannot account for dephasing + + """ + if decoherence: + raise Exception("Decoherence not yet implemented.") + + rhoi = numpy.zeros((self.dim, self.dim), + dtype=REAL) + rhot = ReducedDensityMatrixEvolution(timeaxis=self.TimeAxis, rhoi=rhoi) + for ii in range(self.dim): + for jj in range(self.dim): + rhot.data[:,ii,jj] = self.data[:,ii]*self.data[:,jj] + + return rhot + \ No newline at end of file diff --git a/quantarhei/qm/propagators/oqssvpropagator.py b/quantarhei/qm/propagators/oqssvpropagator.py new file mode 100644 index 0000000..fcc95ac --- /dev/null +++ b/quantarhei/qm/propagators/oqssvpropagator.py @@ -0,0 +1,132 @@ +# -*- coding: utf-8 -*- +import time +import numpy +from numba import njit + +from .oqssvevolution import OQSStateVectorEvolution +from quantarhei.qm import TDRedfieldRateMatrix +import quantarhei as qr + + +class OQSStateVectorPropagator(): + + + def __init__(self, timeaxis=None, current_matrix=None, agg=None, + theory_type="Redfield"): + + + self.TimeAxis = timeaxis + self.Nt = self.TimeAxis.length + self.Odt = self.TimeAxis.data[1]-self.TimeAxis.data[0] + self.dt = self.Odt + self.Nref = 1 + self.KK = current_matrix + + self.theory_type = theory_type + self.agg = agg + + + def setDtRefinement(self, Nref): + """ + The TimeAxis object specifies at what times the propagation + should be stored. We can tell the propagator to use finer + time step for the calculation by setting the refinement. The + refinement is an integer by which the TimeAxis time step should + be devided to get the finer time step. In the code below, we + have dt = 10 in the TimeAxis, but we want to calculate with + dt = 1 + + >>> KK =numpy.array([[0.0, 0.0],[0.0,1.0]], dtype=float) + >>> times = TimeAxis(0,1000,10.0) + >>> pr = OQSStateVectorPropagator(times, KK) + >>> pr.setDtRefinement(10) + + """ + self.Nref = Nref + self.dt = self.Odt/self.Nref + + + def propagate(self, psii, L=4): + + """Short expansion of an exponention to integrate equations of motion + + + Propagation with Hamiltonian only + + + """ + + pr = OQSStateVectorEvolution(self.TimeAxis, psii) + + _CALC(self.KK, psii.data, self.Nt, self.Nref, L, self.dt, pr.data) + + return pr + + + def get_secular_dynamics(self, psii, L=4): + """Here we obtain the secular dynamics for self-consistent methodology + + """ + + pr = OQSStateVectorEvolution(self.TimeAxis, psii) + + if self.theory_type == "Redfield": + + ham = self.agg.get_Hamiltonian() + sbi = self.agg.get_SystemBathInteraction() + + + RR = TDRedfieldRateMatrix(ham, sbi) + + ppi = psii.data*psii.data + + + ppt = numpy.zeros((self.Nt, ham.dim), dtype=qr.REAL) + ppt[0,:] = ppi + + + #RR0 = numpy.zeros((ham.dim, ham.dim, self.Nt)) + #for it in range(self.Nt): + # RR0[:,:,it] = RR.data[it,:,:] + + #t1 = time.time() + _CALC(RR.data, ppi, self.Nt, self.Nref, L, self.dt, ppt) + #_CALC(RR0, ppi, self.Nt, self.Nref, L, self.dt, ppt) + #t2 = time.time() + #print("Done in ", t2-t1, "sec") + pr.data = numpy.sqrt(ppt) + + return pr + + + def get_c0(self, psii, L=4): + """Returns the zero's order solution to the iterative problem + + """ + return self.get_secular_dynamics(psii, L) + + +#@njit(cache=True) +def _CALC(KK, psii, Nt, Nref, L, dt, out): + """Propagation numerics + + """ + + psi1 = psii + psi2 = psii + + indx = 1 + for ii in range(1, Nt): + + for jj in range(Nref): + + for ll in range(1,L+1): + + psi1 = (dt/ll)*numpy.dot(KK[ii, :,:],psi1) + + psi2 += psi1 + psi1 = psi2 + + out[indx,:] = psi2 + indx += 1 + \ No newline at end of file diff --git a/tests/unit/qm/hilbertspace/oqsstatevector_test.py b/tests/unit/qm/hilbertspace/oqsstatevector_test.py index c7e3a0e..7838dcd 100644 --- a/tests/unit/qm/hilbertspace/oqsstatevector_test.py +++ b/tests/unit/qm/hilbertspace/oqsstatevector_test.py @@ -14,6 +14,7 @@ from quantarhei import OQSStateVector from quantarhei import REAL +import quantarhei as qr class TestOQSStateVector(unittest.TestCase): """Tests for the statevector package @@ -62,6 +63,27 @@ def test_creation_from_list(self): self.assertAlmostEqual(psi.puredot(psi), psi.norm()**2) + + def test_evolution(self): + + import numpy + + KK = numpy.array([[-1.0/200.0, 1.0/200.], + [1.0/200.0, -1.0/200.0]], + dtype=REAL) + + timea = qr.TimeAxis(0.0, 1000, 1.0) + + prop = qr.OQSStateVectorPropagator(timea, KK) + + psii = qr.OQSStateVector(data=[numpy.sqrt(0.1), + numpy.sqrt(0.9)]) + + + + psit = prop.propagate(psii) + +