Skip to content

Commit

Permalink
Update to enable varational wavefucntion method
Browse files Browse the repository at this point in the history
  • Loading branch information
@tmancal74
tmancal74 committed Oct 18, 2023
1 parent a6f8355 commit 58099ae
Show file tree
Hide file tree
Showing 2 changed files with 117 additions and 54 deletions.
165 changes: 112 additions & 53 deletions quantarhei/builders/aggregate_base.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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 \
Expand All @@ -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):
Expand All @@ -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
Expand Down Expand Up @@ -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
"""
Expand All @@ -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 + \
Expand All @@ -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
Expand Down Expand Up @@ -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 + \
Expand All @@ -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
Expand Down Expand Up @@ -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
---------------
Expand Down Expand Up @@ -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)

Expand Down
6 changes: 5 additions & 1 deletion quantarhei/qm/hilbertspace/dmoment.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

from .operators import SelfAdjointOperator
from ...core.managers import BasisManaged
from ... import REAL

import numpy
#import scipy
Expand All @@ -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():
Expand Down

0 comments on commit 58099ae

Please sign in to comment.