Some improments and corrections of Foerster theory implementation

examples/symbolic/test_symbolic_fstr.py

@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+from quantarhei.symbolic.cumulant import Uged, Ugde, Uedg, Uegd
+from quantarhei.symbolic.cumulant import gg
+from quantarhei.symbolic.cumulant import CumulantExpr
+from quantarhei.symbolic.abc import a, b, c, d, e, t, T, tau
+
+from sympy import Symbol
+
+""" 
+Cummulant expression for the time-non-local Förster theory
+
+    U_a(t-tau)U_c(tau) W_eq Dagger(U_d(tau))Dagger(U_b(t-tau))
+    
+    =         Dagger(U_d(tau))   Dagger(U_b(t-tau))            U_a(t-tau) U_c(tau) W_eq
+    
+    =[U_g(tau)Dagger(U_d(tau))] [Dagger(U_b(t-tau))U_g(t-tau)]
+     [Dagger(U_g(t-tau))U_a(t-tau)] [U_c(tau)Dagger(U_g(tau))]W_eq
+    
+    
+    = Uged(d,tau) * Uedg(b,t-tau) * Ugde(a,t-tau) * Uegd(c,tau)
+
+"""
+
+
+A = Uged(d,tau)*Uedg(b,t-tau)*Ugde(a,t-tau)*Uegd(c,tau)
+Anorm = Uged(d,tau)*Uegd(c,tau)  
+
+verbatim = True
+
+if verbatim:
+    print(" ")
+    print("Expression to evaluate: ")
+    print(" ")
+    print("    Tr_bath{",A,"W_eq}")
+    print(" ")
+    print("The expression is normalized by:")
+    print(" ")
+    print("    Tr_bath{",Anorm,"W_eq}")
+    print(" ")
+
+
+A = A.rewrite(gg)
+expr = CumulantExpr(A)
+""" use option large=T to evaluate in T --> oo """
+expr = expr.evaluate(large=T) 
+""" use the symetry of lineshape function in the exciton indices """
+#D = CumulantExpr(expr)._leading_index(a)
+#expr = D._getExpr()
+
+A = Anorm.rewrite(gg)
+norm = CumulantExpr(A)
+""" use option large=T to calculate evaluate in T --> oo """
+norm = norm.evaluate(large=T) 
+""" use the symetry of lineshape function in the exciton indices """
+#D = CumulantExpr(expr)._leading_index(a)
+#expr = D._getExpr()
+
+expr = (expr-norm).simplify()
+
+if verbatim:     
+    print("\nCumulant: ")
+    print(" ")
+    print(expr)
+    print(" ")
+
+print("Foerster rate")
+#expr = expr.subs(e,d).simplify()
+
+#a1 = Symbol('a1')
+#b1 = Symbol('b1')
+#expr = expr.subs({d:b1,c:a1,b:a1,a:b1})
+#expr = expr.subs({a1:a,b1:b}).simplify()
+
+print(expr)

quantarhei/builders/molecules.py

@@ -241,33 +241,33 @@ def build(self):
     # I am systematically removing any "naming" in Quantarhei
     #
 
-    # def get_name(self):
-    #     """Returns the name of the molecule
+    def get_name(self):
+        """Returns the name of the molecule
         
-    #     Examples
-    #     --------
+        Examples
+        --------
         
-    #     >>> m = Molecule([0.0, 1.0], name="Jane")
-    #     >>> m.get_name()
-    #     'Jane'
+        >>> m = Molecule([0.0, 1.0], name="Jane")
+        >>> m.get_name()
+        'Jane'
         
-    #     """
-    #     return self.name
+        """
+        return self.name
 
 
-    # def set_name(self, name):
-    #     """Sets the name of the Molecule object
+    def set_name(self, name):
+        """Sets the name of the Molecule object
         
-    #     Examples
-    #     --------
+        Examples
+        --------
         
-    #     >>> m = Molecule([0.0, 1.0])
-    #     >>> m.set_name("Jane")
-    #     >>> print(m.get_name())
-    #     Jane
+        >>> m = Molecule([0.0, 1.0])
+        >>> m.set_name("Jane")
+        >>> print(m.get_name())
+        Jane
         
-    #     """
-    #     self.name = name
+        """
+        self.name = name
 
     def set_electronic_rwa(self, rwa_indices):
         """Sets which electronic states belong to different blocks

quantarhei/builders/opensystem.py

@@ -244,6 +244,7 @@ def get_RelaxationTensor(self, timeaxis,
                 for i in range(ham.dim):
                     dat[i,i] = ham._data[i,i]
                 ham_0 = Hamiltonian(data=dat)
+                ham_0.set_rwa(ham.rwa_indices)
 
             else:
 
@@ -253,11 +254,13 @@ def get_RelaxationTensor(self, timeaxis,
                 # This is done strictly in site basis
                 #
 
-                relaxT = FoersterRelaxationTensor(ham, sbi)
+                relaxT = FoersterRelaxationTensor(ham, sbi, 
+                                                  pure_dephasing=True)
                 dat = numpy.zeros((ham.dim,ham.dim),dtype=REAL)
                 for i in range(ham.dim):
                     dat[i,i] = ham._data[i,i]
                 ham_0 = Hamiltonian(data=dat)
+                ham_0.set_rwa(ham.rwa_indices)
 
             # The Hamiltonian for propagation is the one without
             # resonance coupling
@@ -272,9 +275,12 @@ def get_RelaxationTensor(self, timeaxis,
 
             if time_dependent:
 
+
                 # Time dependent standard Foerster
                 relaxT = NEFoersterRelaxationTensor(ham, sbi, 
                                             as_kernel=as_convolution_kernel)
+
+
                 dat = numpy.zeros((ham.dim,ham.dim),dtype=REAL)
 
                 #for i in range(ham.dim):

quantarhei/qm/hilbertspace/statevector.py

@@ -11,7 +11,7 @@
 
 from ...utils.types import BasisManagedComplexArray
 from ...core.managers import BasisManaged
-from ... import REAL
+from ... import REAL, COMPLEX
 
 class StateVector(BasisManaged):
     """Represents a quantum mechanical state vector
@@ -68,7 +68,7 @@ def __init__(self, dim=None, data=None):
             # check and save dim
             if dim is not None:
                 self.dim = dim
-                self.data = numpy.zeros(dim, dtype=REAL)
+                self.data = numpy.zeros(dim, dtype=COMPLEX)
                 self._initialized = True
 
 

quantarhei/qm/liouvillespace/foerstertensor.py

@@ -15,7 +15,8 @@ class FoersterRelaxationTensor(RelaxationTensor):
     
     
     """
-    def __init__(self, ham, sbi, initialize=True, cutoff_time=None):
+    def __init__(self, ham, sbi, initialize=True, cutoff_time=None, 
+                 pure_dephasing=False):
 
         self._initialize_basis()
 
@@ -29,6 +30,7 @@ def __init__(self, ham, sbi, initialize=True, cutoff_time=None):
         self._is_initialized = False
         self._has_cutoff_time = False
         self.as_operators = False
+        self.pure_dephasing = pure_dephasing
 
         if cutoff_time is not None:
             self.cutoff_time = cutoff_time
@@ -83,3 +85,31 @@ def initialize(self):
             # calculate dephasing rates and depopulation rates
             #
             self.updateStructure()
+
+            if self.pure_dephasing:
+                # additional pure dephasing 
+                self.add_dephasing()
+
+
+    def add_dephasing(self):
+
+
+        # line shape function derivatives
+        sbi = self.SystemBathInteraction
+        Na = self.dim
+        Nt = sbi.TimeAxis.length
+
+        ht = numpy.zeros((Na, Nt),
+                         dtype=numpy.complex64)
+
+        sbi.CC.create_one_integral()
+
+        for ii in range(1, Na):
+           ht[ii,:] = sbi.CC.get_hoft(ii-1,ii-1)
+
+
+        for aa in range(self.dim):
+            for bb in range(self.dim):
+                if aa != bb:
+
+                    self.data[aa,bb,aa,bb] -= (ht[aa,Nt-1]+ht[bb,Nt-1])