From 58099ae7f28ad3b82f8de73294fe588e89e8ed35 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Man=C4=8Dal?= Date: Wed, 18 Oct 2023 14:00:56 +0200 Subject: [PATCH] Update to enable varational wavefucntion method --- quantarhei/builders/aggregate_base.py | 165 +++++++++++++++++--------- quantarhei/qm/hilbertspace/dmoment.py | 6 +- 2 files changed, 117 insertions(+), 54 deletions(-) diff --git a/quantarhei/builders/aggregate_base.py b/quantarhei/builders/aggregate_base.py index e834191..9d9aea1 100644 --- a/quantarhei/builders/aggregate_base.py +++ b/quantarhei/builders/aggregate_base.py @@ -1327,7 +1327,7 @@ def __str__(self): def build(self, mult=1, sbi_for_higher_ex=False, vibgen_approx=None, Nvib=None, vibenergy_cutoff=None, - fem_full=False): + fem_full=False, el_blocks=False): """Builds aggregate properties Calculates Hamiltonian and transition dipole moment matrices and @@ -1423,6 +1423,15 @@ def build(self, mult=1, sbi_for_higher_ex=False, vibenergy_cutoff=vibenergy_cutoff): self.all_states.append((a, s1)) + if el_blocks: + self.electronic_blocks = dict() + for sig in self.elsigs: + self.electronic_blocks[sig] = "ciao" + + + print("AHOJ !!!!!!!!!!!!!!!!!") + print(self.electronic_blocks) + # Set up Hamiltonian and Transition dipole moment matrices for a, s1 in self.all_states: #self.allstates(mult=self.mult, @@ -1765,6 +1774,12 @@ def convert_to_ground_vibbasis(self, operator, Nt=None): nop = Hamiltonian(dim=self.Ntot) nop1 = Hamiltonian(dim=self.Nel) + + + elif isinstance(operator, TransitionDipoleMoment): + + nop = TransitionDipoleMoment(dim=self.Ntot) + n_indices = 3 elif isinstance(operator, ReducedDensityMatrixEvolution) or \ @@ -1786,46 +1801,68 @@ def convert_to_ground_vibbasis(self, operator, Nt=None): else: raise Exception("Operation not implemented for this type: "+ operator.__class__.__name__) + + # 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") + + # + # 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) 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") + # # 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): @@ -1840,26 +1877,43 @@ def convert_to_ground_vibbasis(self, operator, Nt=None): 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] + self.FCf[i_ng, j_n]*\ + operator._data[j_n, j_m]*\ + self.FCf[j_m, i_mg] + return nop #, nop1 + + elif n_indices == 3: + + # do the conversion + + #print("TRANSFORMING") + # loop over 3D + for a in range(3): + # 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,a] = 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, a]*\ + self.FCf[j_m, i_mg] + + return nop + else: raise Exception("Incompatible operator") + def trace_converted(self, operator, Nt=None): n_indices = 2 @@ -2562,7 +2616,7 @@ def _thermal_population(self, temp=0.0, subtract=None, def _impulsive_population(self, relaxation_theory_limit="weak_coupling", - temperature=0.0): + temperature=0.0, DD=None): """Impulsive excitation of the density matrix from ground state """ @@ -2572,7 +2626,8 @@ def _impulsive_population(self, relaxation_theory_limit="weak_coupling", temperature=temperature) rho0 = rho.data - DD = self.TrDMOp.data + if DD is None: + DD = self.TrDMOp.data # abs value of the transition dipole moment dabs = numpy.sqrt(DD[:,:,0]**2 + \ @@ -2583,7 +2638,7 @@ def _impulsive_population(self, relaxation_theory_limit="weak_coupling", return rho0 - def get_StateVector(self, condition_type=None): + def get_StateVector(self, condition_type=None, DD=None): """Returns state vector accordoing to specified conditions Parameters @@ -2614,7 +2669,8 @@ def get_StateVector(self, condition_type=None): elif condition_type == "impulsive_excitation": - DD = self.TrDMOp.data + if DD is None: + DD = self.TrDMOp.data # abs value of the transition dipole moment dabs = numpy.sqrt(DD[:,:,0]**2 + \ @@ -2637,7 +2693,7 @@ def get_StateVector(self, condition_type=None): def get_DensityMatrix(self, condition_type=None, relaxation_theory_limit="weak_coupling", temperature=None, - relaxation_hamiltonian=None): + relaxation_hamiltonian=None, DD=None): """Returns density matrix according to specified condition Returs density matrix to be used e.g. as initial condition for @@ -2665,6 +2721,9 @@ def get_DensityMatrix(self, condition_type=None, Hamiltonian according to which we form thermal equilibrium. In case of `strong_coupling`, no reorganization energies are subtracted - we assume that the supplied energies are already void of them. + + DD : real array + Submit your own transition dipole moment matrix (for x, y and z) Condition types --------------- @@ -2705,7 +2764,7 @@ def get_DensityMatrix(self, condition_type=None, elif condition_type == "impulsive_excitation": rho0 = self._impulsive_population( relaxation_theory_limit=relaxation_theory_limit, - temperature=temperature) + temperature=temperature, DD=DD) self.rho0 = rho0 return DensityMatrix(data=self.rho0) diff --git a/quantarhei/qm/hilbertspace/dmoment.py b/quantarhei/qm/hilbertspace/dmoment.py index 0aac518..70cdf4c 100644 --- a/quantarhei/qm/hilbertspace/dmoment.py +++ b/quantarhei/qm/hilbertspace/dmoment.py @@ -2,6 +2,7 @@ from .operators import SelfAdjointOperator from ...core.managers import BasisManaged +from ... import REAL import numpy #import scipy @@ -19,7 +20,10 @@ def __init__(self, dim=None, data=None): if cb != 0: self.manager.register_with_basis(cb,self) - self._data = data + if data is None: + self._data = numpy.zeros((dim,dim,3), dtype=REAL) + else: + self._data = data self.dim = self._data.shape[0] if not self.check_selfadjoint():