Skip to content

Commit

Permalink
Various updates
Browse files Browse the repository at this point in the history 
First attempts on non-linear wavefunction equations
  • Loading branch information
@tmancal74
tmancal74 committed Oct 16, 2023
1 parent 9880031 commit a6f8355
Show file tree
Hide file tree
Showing 4 changed files with 395 additions and 8 deletions.
238 changes: 238 additions & 0 deletions quantarhei/builders/aggregate_base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1737,6 +1737,188 @@ def rebuild(self, mult=1, sbi_for_higher_ex=False,
#
###########################################################################

def convert_to_ground_vibbasis(self, operator, Nt=None):
"""Converts an operator to a ground state vibrational basis repre
Default representation in Quantarhei is that with a specific shifted
vibrational basis in each electronic state. Here we convert to a
representation where there is a single basis used for all vibrational
states regardless of elecrtronic state
The conversion MUST be done in site basis. Only in site basis
we can distinguish the vibrational states properly
"""
n_indices = 2
evolution = False
whole = False
if operator.dim == self.Ntot:

if isinstance(operator, ReducedDensityMatrix) or \
isinstance(operator, DensityMatrix):

nop = ReducedDensityMatrix(dim=self.Ntot)


elif isinstance(operator, Hamiltonian):

nop = Hamiltonian(dim=self.Ntot)
nop1 = Hamiltonian(dim=self.Nel)


elif isinstance(operator, ReducedDensityMatrixEvolution) or \
isinstance(operator, DensityMatrixEvolution):

if Nt is not None:
nop = ReducedDensityMatrix(dim=self.Ntot)
evolution = True
whole = False
else:
nop = ReducedDensityMatrixEvolution(operator.TimeAxis)
rhoi = ReducedDensityMatrix(dim=self.Ntotl)
# we set zero initial condition, because this initialized
# the data storage
nop.set_initial_condition(rhoi)
evolution = True
whole = True

else:
raise Exception("Operation not implemented for this type: "+
operator.__class__.__name__)

if n_indices == 2:

# convert to representation by ground-state oscillator

# FIXME: This limitation might not be necessary
# in the ground states of all monomers, there must be the same
# or greater number of levels than in the excited state

# over all monomers
for k in range(self.nmono):
mono = self.monomers[k]
# over all modes
n_mod = mono.get_number_of_modes()
for i in range(n_mod):
mod = mono.get_Mode(i)
n_g = mod.get_nmax(0)
# over all excited states
# FIXME: this should be mono.Nel as in Aggregate
for j in range(mono.nel):
if (j > 0):
n_e = mod.get_nmax(j)
if n_e > n_g:
raise Exception("Number of levels"+
" in the excited state of a molecule has to be \n"+
"the same or smaller than in the ground state")


# do the conversion

#
# ground state vibrational states
#
stgs = []
for i_g in self.vibindices[0]:
vs_g = self.vibsigs[i_g]
stg = self.get_VibronicState(vs_g[0],
vs_g[1])
stgs.append(stg)

#print("TRANSFORMING")
# loop over electronic states n, m
for n in range(self.Nel):
i_ng = -1
for i_n in self.vibindices[n]:
i_ng += 1
for m in range(self.Nel):
i_mg = -1
for i_m in self.vibindices[m]:
i_mg += 1
nop._data[i_n,i_m] = 0.0
for j_n in self.vibindices[n]:
for j_m in self.vibindices[m]:
nop._data[i_n,i_m] += \
self.FCf[i_ng, j_n]*operator._data[j_n, j_m]*self.FCf[j_m, i_mg]
#if n == 1 and m == 1:
# print(self.FCf[i_n, j_n],self.FCf[j_m, i_m],i_n,j_n, j_m, i_m)

#for n in range(self.Nel):
# i_ng = -1
# for i_n in self.vibindices[n]:
# i_ng += 1
# for m in range(self.Nel):
# i_mg = -1
# for i_m in self.vibindices[m]:
# i_mg += 1
# if i_ng == i_mg:
# nop1._data[n,m] += \
# nop._data[i_n, i_m]
return nop #, nop1

else:
raise Exception("Incompatible operator")

def trace_converted(self, operator, Nt=None):

n_indices = 2
evolution = False
whole = False

if operator.dim == self.Ntot:

if isinstance(operator, ReducedDensityMatrix) or \
isinstance(operator, DensityMatrix):

nop = ReducedDensityMatrix(dim=self.Nel)


elif isinstance(operator, Hamiltonian):

nop = Hamiltonian(dim=self.Nel)


elif isinstance(operator, ReducedDensityMatrixEvolution) or \
isinstance(operator, DensityMatrixEvolution):

if Nt is not None:
nop = ReducedDensityMatrix(dim=self.Nel)
evolution = True
whole = False
else:
nop = ReducedDensityMatrixEvolution(operator.TimeAxis)
rhoi = ReducedDensityMatrix(dim=self.Nel)
# we set zero initial condition, because this initialized
# the data storage
nop.set_initial_condition(rhoi)
evolution = True
whole = True

else:
raise Exception("Operation not implemented for this type: "+
operator.__class__.__name__)

if n_indices == 2:


for n in range(self.Nel):
i_ng = -1
for i_n in self.vibindices[n]:
i_ng += 1
for m in range(self.Nel):
i_mg = -1
for i_m in self.vibindices[m]:
i_mg += 1
if i_ng == i_mg:
nop._data[n,m] += \
operator._data[i_n, i_m]
return nop

else:
raise Exception("Incompatible operator")


def trace_over_vibrations(self, operator, Nt=None):
"""Average an operator over vibrational degrees of freedom
Expand All @@ -1756,6 +1938,11 @@ def trace_over_vibrations(self, operator, Nt=None):
nop = ReducedDensityMatrix(dim=self.Nel)


elif isinstance(operator, Hamiltonian):

nop = Hamiltonian(dim=self.Nel)


elif isinstance(operator, ReducedDensityMatrixEvolution) or \
isinstance(operator, DensityMatrixEvolution):

Expand Down Expand Up @@ -2396,6 +2583,57 @@ def _impulsive_population(self, relaxation_theory_limit="weak_coupling",
return rho0


def get_StateVector(self, condition_type=None):
"""Returns state vector accordoing to specified conditions
Parameters
----------
condition_type : str
Type of the initial condition. If None, the property sv0, which
was presumably calculated in the past, is returned.
Condition types
---------------
impulsive_excitation
Excitation by ultrabroad laser pulse
"""
# aggregate must be built before we call this method
if not self._built:
raise Exception("Aggregate must be built before"
+" get_StateVector can be invoked.")

# if no condition is specified, it is understood that we return
# internal sv0, which was calculated sometime in the past
if condition_type is None:
return StateVector(data=self.sv0)


elif condition_type == "impulsive_excitation":

DD = self.TrDMOp.data

# abs value of the transition dipole moment
dabs = numpy.sqrt(DD[:,:,0]**2 + \
DD[:,:,1]**2 + DD[:,:,2]**2)

sv0 = numpy.zeros(self.Ntot, dtype=COMPLEX)

# zero temperature
sv0[0] = 1.0

self.sv0 = numpy.dot(dabs,sv0)

return StateVector(data=self.sv0)

else:
raise Exception("Unknown condition type")



def get_DensityMatrix(self, condition_type=None,
relaxation_theory_limit="weak_coupling",
temperature=None,
Expand Down
24 changes: 19 additions & 5 deletions quantarhei/qm/propagators/rdmpropagator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ class ReducedDensityMatrixPropagator(MatrixData, Saveable):
"""

def __init__(self, timeaxis=None, Ham=None, RTensor=None, Iterm=None,
Efield=None, Trdip=None, PDeph=None):
Efield=None, Trdip=None, PDeph=None, NonHerm=None):
"""
Creates a Reduced Density Matrix propagator which can propagate
Expand Down Expand Up @@ -94,6 +94,7 @@ def __init__(self, timeaxis=None, Ham=None, RTensor=None, Iterm=None,
self.has_Iterm = False
self.has_RWA = False
self.has_EField = False
self.has_NonHerm = False

if not ((timeaxis is None) and (Ham is None)):

Expand Down Expand Up @@ -166,6 +167,9 @@ def __init__(self, timeaxis=None, Ham=None, RTensor=None, Iterm=None,
else:
raise Exception("RelaxationTensor has to be set first.")

if NonHerm is not None:
self.has_NonHerm = True
self.NonHerm = NonHerm

self.Odt = self.TimeAxis.data[1]-self.TimeAxis.data[0]
self.dt = self.Odt
Expand Down Expand Up @@ -495,7 +499,8 @@ def __propagate_short_exp(self, rhoi, L=4):

for ll in range(1,L+1):

rho1 = -_COM(HH, ll, self.dt, rho1)
rho1 = -_COM(HH, ll, self.dt, rho1,
has_NonHerm=self.has_NonHerm)

rho2 = rho2 + rho1
rho1 = rho2
Expand Down Expand Up @@ -550,7 +555,12 @@ def _INIT_RWA(self):
if self.Hamiltonian.has_rwa:
HH = self.Hamiltonian.get_RWA_data()
else:
HH = self.Hamiltonian.data
HH = self.Hamiltonian.data


if self.has_NonHerm:
HH = HH + 1j*self.NonHerm

return HH


Expand Down Expand Up @@ -1603,7 +1613,7 @@ def _OTI(rhoY, Km, Kd, Lm, Ld, ll, dt, rho1):
)


def _COM(HH, ll, dt, rho1):
def _COM(HH, ll, dt, rho1, has_NonHerm=False):
"""Commutator with the a given operator (e.g. Hamiltonian)
Parameters
Expand All @@ -1626,7 +1636,11 @@ def _COM(HH, ll, dt, rho1):
"""
ret = (1j*dt/ll)*(numpy.dot(HH,rho1) - numpy.dot(rho1,HH))
if has_NonHerm:
H2 = numpy.conj(numpy.transpose(HH))
else:
H2 = HH
ret = (1j*dt/ll)*(numpy.dot(HH,rho1) - numpy.dot(rho1,H2))
return ret


Expand Down
Loading

0 comments on commit a6f8355

Please sign in to comment.