From 102e47e96894d4fe41afc8d429649bbad2bccd1a Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 23 Jun 2023 15:52:24 +0200 Subject: [PATCH 01/25] Introducing attributes|properties 'classical' and 'bipartite' to classify system Hamiltonians --- wave_train/hamilton/chain.py | 8 ++++++-- wave_train/hamilton/coupled.py | 2 ++ wave_train/hamilton/exciton.py | 2 ++ wave_train/hamilton/phonon.py | 2 ++ 4 files changed, 12 insertions(+), 2 deletions(-) diff --git a/wave_train/hamilton/chain.py b/wave_train/hamilton/chain.py index ce79120..6fabd5d 100644 --- a/wave_train/hamilton/chain.py +++ b/wave_train/hamilton/chain.py @@ -22,7 +22,10 @@ def __init__(self, n_site, periodic, homogen): Homogeneous chain/ring really_homogen: boolean Really a homogeneous chain/ring - + classical: boolean + Whether a system has a well-defined classical limit + bipartite: boolean + Whether a system is bipartite, i.e. comprising two classes of particles. """ if n_site < 2: @@ -32,8 +35,9 @@ def __init__(self, n_site, periodic, homogen): self.n_site = n_site self.periodic = periodic self.homogen = homogen - self.really_homogen = self.homogen + self.classical = [] + self.bipartite = [] def __str__(self): diff --git a/wave_train/hamilton/coupled.py b/wave_train/hamilton/coupled.py index 3323735..3d2c581 100644 --- a/wave_train/hamilton/coupled.py +++ b/wave_train/hamilton/coupled.py @@ -76,6 +76,8 @@ def __init__(self, n_site, periodic, homogen, # Construct object of parent class Chain.__init__(self, n_site, periodic, homogen) + self.classical = True + self.bipartite = True # Construct objects of phonon and exciton class self.ex = Exciton(n_site, periodic, homogen, alpha, beta, eta) diff --git a/wave_train/hamilton/exciton.py b/wave_train/hamilton/exciton.py index e520f8e..7e91b7b 100644 --- a/wave_train/hamilton/exciton.py +++ b/wave_train/hamilton/exciton.py @@ -35,6 +35,8 @@ def __init__(self, n_site, periodic, homogen, alpha=1, beta=1, eta=0): # Initialize object of parent class Chain.__init__(self, n_site, periodic, homogen) + self.classical = False + self.bipartite = False # Parameters of a homogeneous chain/ring are converted from scalar to vector if homogen: diff --git a/wave_train/hamilton/phonon.py b/wave_train/hamilton/phonon.py index 8eb9444..3c8f01b 100644 --- a/wave_train/hamilton/phonon.py +++ b/wave_train/hamilton/phonon.py @@ -42,6 +42,8 @@ def __init__(self, n_site, periodic, homogen, mass=1, nu=1, omg=1): # Construct object of parent class Chain.__init__(self, n_site, periodic, homogen) + self.classical = True + self.bipartite = False # Parameters of a homogeneous chain/ring are converted from scalar to vector if homogen: From c15275a6c2fcf2f38d76377b0f58738d0fe09357 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 23 Jun 2023 16:54:37 +0200 Subject: [PATCH 02/25] Start doing some experiments beyond Exciton, Phonon, Coupled --- test_scripts/XYZ/tise_1.py | 75 +++++++++ wave_train/hamilton/xyz.py | 306 +++++++++++++++++++++++++++++++++++++ 2 files changed, 381 insertions(+) create mode 100644 test_scripts/XYZ/tise_1.py create mode 100644 wave_train/hamilton/xyz.py diff --git a/test_scripts/XYZ/tise_1.py b/test_scripts/XYZ/tise_1.py new file mode 100644 index 0000000..4516d0e --- /dev/null +++ b/test_scripts/XYZ/tise_1.py @@ -0,0 +1,75 @@ +from ...wave_train.hamilton.exciton import Exciton +from ...wave_train.dynamics.tise import TISE +from ...wave_train.io.logging import TeeLogger +from os.path import basename, splitext + +def xyz_tise(batch_mode): + # Detect name of this script file (without extension) + base_name = basename(__file__) + my_file = splitext(base_name)[0] + + # logging instance: will be initialized with + # class for logging to both console and logfile + logger = None + if not batch_mode: + logger = TeeLogger(log_file=my_file + ".log") + + # Set up the excitonic Hamiltonian for a chain + hamilton = Exciton( + n_site=6, # number of sites + periodic=True, # periodic boundary conditions + homogen=True, # homogeneous chain/ring + alpha=0.1, # excitonic site energy + beta=-0.01, # coupling strength (NN) + eta=0.0 # constant energy offset + ) + + # Set up TT representation of the Hamiltonian + hamilton.get_TT( + n_basis=2, # size of electronic basis set + qtt=False # using quantized TT format + ) + + # Set up TISE solver + # n_levels = n_site+1 # ground and singly excited + # n_levels = n_site*(n_site+1)//2+1 # ground and singly and doubly excited + dynamics = TISE( + hamilton=hamilton, # choice of Hamiltonian, see above + n_levels=16, # number of energy levels to be calculated + solver='als', # choice of eigensolver for the full system + eigen='eig', # choice of eigensolver for the micro systems + ranks=15, # rank of initial guess for ALS + repeats=20, # number of sweeps in eigensolver scheme + conv_eps=1e-8, # threshold for detecting convergence of the eigenvalue + e_est=0.0, # estimation: eigenvalues closest to this number + e_min=0.05, # lower end of energy plot axis (if exact energies not available!) + e_max=0.15, # upper end of energy plot axis (if exact energies not available!) + save_file=my_file+'.pic', # if not None, generated data will be saved to this file + load_file=None, # if not None, reference data will be loaded from this file + compare=None # type of comparison with reference data + ) + + # Batch mode + if batch_mode: + dynamics.solve() # Solve TISE *without* visualization + + # Interactive mode: Setup animated visualization + else: + from wave_train.graphics.factory import VisualTISE + graphics = VisualTISE( + dynamics=dynamics, # choice of dynamics (EoM), see above + plot_type='QuantNumbers', # select your favorite plot type + plot_expect=True, # toggle plotting of expectation values + figure_pos=(100, 50), # specifying position (x,y) of upper left of figure [in pixels] + figure_size=(1050, 450), # specifying size (w,h) of figure [in pixels] + image_file=my_file+'.png', # if not None, image (last frame) will be written to this file + movie_file=my_file+'.mp4', # if not None, animation will be written to this file + snapshots=False, # each snapshot is saved + frame_rate=1, # frames per second in mp4 animation file + plot_style={} # additional plot style information + ).create() + graphics.solve() # Solve TISE *with* visualization + + +if __name__ == '__main__': + xyz_tise(batch_mode=False) diff --git a/wave_train/hamilton/xyz.py b/wave_train/hamilton/xyz.py new file mode 100644 index 0000000..37f2d2e --- /dev/null +++ b/wave_train/hamilton/xyz.py @@ -0,0 +1,306 @@ +import sys +import numpy as np +from wave_train.hamilton.chain import Chain + + +class XYZ(Chain): + + """ + Exciton dynamics in one dimension + --------------------------------- + + Linear chain of electronic sites with periodic boundaries + Model consists of site energies (alpha) and nearest neighbor coupling (beta) + + Parent class is chain. + """ + + def __init__(self, n_site, periodic, homogen, alpha=1, beta=1, eta=0): + """ + Parameters: + ----------- + n_site: int + Length of the chain + periodic: boolean + Periodic boundary conditions + homogen: boolean + Homogeneous chain/ring + alpha: float + excitonic site energy + beta: float + coupling strength between nearest neighbours + eta: float + constant energy offset + """ + + # Initialize object of parent class + Chain.__init__(self, n_site, periodic, homogen) + self.classical = False + self.bipartite = False + + # Parameters of a homogeneous chain/ring are converted from scalar to vector + if homogen: + if isinstance(alpha,list): + sys.exit("Wrong input of alpha parameter: should be a scalar") + self.alpha = np.repeat(alpha, n_site) + + if isinstance(beta,list): + sys.exit("Wrong input of beta parameters: should be a scalar") + if periodic: + self.beta = np.repeat(beta, n_site) + else: + self.beta = np.repeat(beta, n_site-1) + + # Parameters of an inhomogeneous chain/ring should be given as vectors + else: + if len(alpha) != n_site: + sys.exit("Inconsistent length of vector of alpha parameters") + self.alpha = alpha + + if periodic: + if len(beta) != n_site: + sys.exit("Inconsistent length of vector of beta parameters") + else: + if len(beta) != n_site - 1: + sys.exit("Inconsistent length of vector of beta parameters") + self.beta = beta + + # Constant energy offset + self.eta = eta + + # Print output from superclass and from this class + print(super().__str__()) + print(self) + + # Useful for savemat|pickle output: Name and classification of this class + self.name = self.__class__.__name__ + self.bipartite = False # Single class of (quasi-) particles + self.classical = False # Classical approximation not supported + + def __str__(self): + if self.homogen: + alpha = self.alpha[0] + beta = self.beta[0] + else: + alpha = self.alpha + beta = self.beta + eta = self.eta + exciton_str = """ +------------------------ +Hamiltonian for EXCITONS +------------------------ + +Excitonic site energy (alpha) = {} +NN coupling strength (beta) = {} +Constant energy offset (eta) = {} + """.format(alpha, beta, eta) + + return exciton_str + + def get_2Q(self, n_dim=2): + """" + Second quantization: + Formulate Hamiltonian in terms of creation/annihiliation operators + """ + + # Call corresponding method from superclass Chain + super().get_2Q(n_dim) + + # Set up ground and first excited state + self.ground = np.zeros(n_dim) + self.ground[0] = 1 + self.excited = np.zeros(n_dim) + self.excited[1] = 1 + + def get_SLIM(self, n_dim=2): + """ + Set up SLIM matrices and operators of second quantization, + mainly for use in tensor train representation + + Parameters: + ----------- + n_dim: int + dimension of truncated(!) Hilbert space + + Returns + ------- + S: ndarray or list of ndarrays + single-site components of SLIM decomposition + L: ndarray or list of ndarrays + left two-site components of SLIM decomposition + I: ndarray or list of ndarrays + identity components of SLIM decomposition + M: ndarray or list of ndarrays + right two-site components of SLIM decomposition + + Reference: + Gelss/Klus/Matera/Schuette : Nearest-neighbor interaction systems in the tensor-train format + https://doi.org/10.1016/j.jcp.2017.04.007 + """ + + # Size of basis set + if n_dim != 2: + sys.exit(""" +--------------------------------- +Value of n_dim must be set to two +--------------------------------- + """) + self.n_dim = n_dim + + # Using second quantization + self.get_2Q (self.n_dim) + + # Set up S,L,I,M operators + S = self.qu_numbr + L_1 = self.raising + L_2 = self.lowering + I = self.identity + M_1 = self.lowering + M_2 = self.raising + + if self.homogen: + + # define ndarrays + self.S = self.alpha[0] * S + self.eta * I / self.n_site + self.L = self.beta[0] * np.stack((L_1, L_2), axis=-1) + self.I = I + self.M = np.stack((M_1, M_2)) + + else: + + # define lists of ndarrays + self.S = [None] * self.n_site + self.L = [None] * self.n_site + self.I = [None] * self.n_site + self.M = [None] * self.n_site + + # insert SLIM components + for i in range(self.n_site): # all sites (S) + self.S[i] = self.alpha[i] * S + self.eta * I / self.n_site + + for i in range(self.n_site - 1): # all sites but last (L) + self.L[i] = self.beta[i] * np.stack((L_1, L_2), axis=-1) + # self.L[i][:, :, 0] = self.beta[i] * L_1 + # self.L[i][:, :, 1] = self.beta[i] * L_2 + + for i in range(self.n_site): # all sites (I) + self.I[i] = I + + for i in range(1, self.n_site): # all sites but first (M) + self.M[i] = np.stack((M_1, M_2)) + # self.M[i][0, :, :] = M_1 + # self.M[i][1, :, :] = M_2 + + if self.periodic: # coupling last site (L) to first site (M) + self.L[self.n_site - 1] = self.beta[self.n_site - 1] * np.stack((L_1, L_2), axis=-1) + self.M[0] = np.stack((M_1, M_2)) + + print(""" +------------------------- +Setting up SLIM operators +------------------------- + +Size of excitonic basis set = {} + """.format(self.n_dim)) + + if self.homogen: + print('shape of S = ', self.S.shape) + print('shape of L = ', self.L.shape) + print('shape of I = ', self.I.shape) + print('shape of M = ', self.M.shape) + + def get_exact (self, n_levels): + """Compute analytic (exact) phononic energy levels for homogeneous systems, + i.e. analytic solutions of the underlying TISE. + So far, only up to two quanta of excitation + + Parameters: + ----------- + n_levels: int + number of exact energies to be computed + + Return: + ------ + energies: list of floats + list of exact energy levels + + Reference: + ---------- + Wikipedia article : Hueckel method + https://en.wikipedia.org/wiki/H%C3%BCckel_method + Z Hu, G. S. Engels, S. Kais : J. Chem. Phys. 148, 204307 (2018) + https://doi.org/10.1063/1.5026116 + """ + + # This method works only for homogeneous chains|rings + if not self.homogen: + return + + # Get necessary parameters + n_site = self.n_site + alpha = self.alpha[0] + beta = self.beta[0] + eta = self.eta + + # index j is an even integer running from 0 to 2N-2 + j = np.linspace(0, 2*n_site-2, n_site) + + # ground state: zero quanta of excitation + energies = eta * np.ones(1) + n_excite = 0 + + if self.periodic: # Cyclic systems (rings) + + # states with one quantum of excitation + if n_levels>1: + n_excite = 1 + ak = np.pi * j / n_site # wavenumber for lattice constant a=1 + ec = alpha + 2 * beta * np.cos(ak) # component energy + energies = np.append(energies, ec + eta) + + # states with two quanta of excitation + if n_levels>1+n_site: + n_excite = 2 + ak = np.pi * (j+1) / n_site # wavenumber for lattice constant a=1 + ec = alpha + 2 * beta * np.cos(ak) # component energy + ee,ff = np.meshgrid(ec,ec) # coordinate matrices from coordinate vectors + em = ee + ff # summing up two component energies + for i in range(1,n_site): # loop over all superdiagonals + energies = np.append(energies, np.diag(em,i)+ eta ) + + # states with more than two quanta not (yet?) implemented + if n_levels>1+sum(range(n_site+1)): + print ('Wrong choice of num_levels : ' + str(n_levels)) + sys.exit(1) + + else: # Linear systems (chains) + + # states with one quantum of excitation + if n_levels>1: + n_excite = 1 + ak = np.pi * (j/2+1) / (n_site+1) + ec = alpha + 2 * beta * np.cos(ak) + energies = np.append(energies, ec) + + # states with more than one quantum not (yet?) implemented + if n_levels>1+n_site: + print ('Wrong choice of num_levels : ' + str(n_levels)) + sys.exit(1) + + # sorting the energies in ascending order + energies = np.sort(energies) + + # truncate unnecessray energy levels + energies = energies[:n_levels] + + print(""" +------------------------------------------------ +Calculating excitonic energy levels analytically +------------------------------------------------ + +Number of energy levels desired : {} +Maximum number of excitations : {} + """.format(len(energies), n_excite)) + + return energies + From 6f544deb4e1a1c6d726eeef5e3882d4682ecca9a Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 23 Jun 2023 16:58:55 +0200 Subject: [PATCH 03/25] Not using "from" in import statements any more --- test_scripts/XYZ/tise_1.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test_scripts/XYZ/tise_1.py b/test_scripts/XYZ/tise_1.py index 4516d0e..f63ff72 100644 --- a/test_scripts/XYZ/tise_1.py +++ b/test_scripts/XYZ/tise_1.py @@ -1,6 +1,6 @@ -from ...wave_train.hamilton.exciton import Exciton -from ...wave_train.dynamics.tise import TISE -from ...wave_train.io.logging import TeeLogger +import wave_train.hamilton.exciton +import wave_train.dynamics.tise +import wave_train.io.logging from os.path import basename, splitext def xyz_tise(batch_mode): @@ -12,10 +12,10 @@ def xyz_tise(batch_mode): # class for logging to both console and logfile logger = None if not batch_mode: - logger = TeeLogger(log_file=my_file + ".log") + logger = wave_train.io.logging.TeeLogger(log_file=my_file + ".log") # Set up the excitonic Hamiltonian for a chain - hamilton = Exciton( + hamilton = wave_train.hamilton.exciton.Exciton( n_site=6, # number of sites periodic=True, # periodic boundary conditions homogen=True, # homogeneous chain/ring @@ -33,7 +33,7 @@ def xyz_tise(batch_mode): # Set up TISE solver # n_levels = n_site+1 # ground and singly excited # n_levels = n_site*(n_site+1)//2+1 # ground and singly and doubly excited - dynamics = TISE( + dynamics = wave_train.dynamics.tise.TISE( hamilton=hamilton, # choice of Hamiltonian, see above n_levels=16, # number of energy levels to be calculated solver='als', # choice of eigensolver for the full system From 564e3e453a585152b756031da47ea8adf964bb66 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 23 Jun 2023 17:03:40 +0200 Subject: [PATCH 04/25] Trying to import xyz class for "fake" Hamiltonian --- test_scripts/XYZ/tise_1.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test_scripts/XYZ/tise_1.py b/test_scripts/XYZ/tise_1.py index f63ff72..30486c8 100644 --- a/test_scripts/XYZ/tise_1.py +++ b/test_scripts/XYZ/tise_1.py @@ -1,4 +1,4 @@ -import wave_train.hamilton.exciton +import wave_train.hamilton.xyz import wave_train.dynamics.tise import wave_train.io.logging from os.path import basename, splitext @@ -15,7 +15,7 @@ def xyz_tise(batch_mode): logger = wave_train.io.logging.TeeLogger(log_file=my_file + ".log") # Set up the excitonic Hamiltonian for a chain - hamilton = wave_train.hamilton.exciton.Exciton( + hamilton = wave_train.hamilton.exciton.XYZ( n_site=6, # number of sites periodic=True, # periodic boundary conditions homogen=True, # homogeneous chain/ring From c44bd79ff7e1838d05648fe46d9f4a35c4420ec9 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Sat, 24 Jun 2023 11:07:18 +0200 Subject: [PATCH 05/25] Undo some of yesterday's changes because attributes "classical" and "bipartite" already existed --- wave_train/hamilton/chain.py | 3 ++- wave_train/hamilton/exciton.py | 2 -- wave_train/hamilton/phonon.py | 2 -- wave_train/hamilton/xyz.py | 4 +++- 4 files changed, 5 insertions(+), 6 deletions(-) diff --git a/wave_train/hamilton/chain.py b/wave_train/hamilton/chain.py index 6fabd5d..56c8813 100644 --- a/wave_train/hamilton/chain.py +++ b/wave_train/hamilton/chain.py @@ -25,7 +25,8 @@ def __init__(self, n_site, periodic, homogen): classical: boolean Whether a system has a well-defined classical limit bipartite: boolean - Whether a system is bipartite, i.e. comprising two classes of particles. + Whether a system is bipartite, i.e. comprising + two classes of particles, e.g. fast and slow. """ if n_site < 2: diff --git a/wave_train/hamilton/exciton.py b/wave_train/hamilton/exciton.py index 7e91b7b..e520f8e 100644 --- a/wave_train/hamilton/exciton.py +++ b/wave_train/hamilton/exciton.py @@ -35,8 +35,6 @@ def __init__(self, n_site, periodic, homogen, alpha=1, beta=1, eta=0): # Initialize object of parent class Chain.__init__(self, n_site, periodic, homogen) - self.classical = False - self.bipartite = False # Parameters of a homogeneous chain/ring are converted from scalar to vector if homogen: diff --git a/wave_train/hamilton/phonon.py b/wave_train/hamilton/phonon.py index 3c8f01b..8eb9444 100644 --- a/wave_train/hamilton/phonon.py +++ b/wave_train/hamilton/phonon.py @@ -42,8 +42,6 @@ def __init__(self, n_site, periodic, homogen, mass=1, nu=1, omg=1): # Construct object of parent class Chain.__init__(self, n_site, periodic, homogen) - self.classical = True - self.bipartite = False # Parameters of a homogeneous chain/ring are converted from scalar to vector if homogen: diff --git a/wave_train/hamilton/xyz.py b/wave_train/hamilton/xyz.py index 37f2d2e..6dd8dec 100644 --- a/wave_train/hamilton/xyz.py +++ b/wave_train/hamilton/xyz.py @@ -8,7 +8,9 @@ class XYZ(Chain): """ Exciton dynamics in one dimension --------------------------------- - + + This is a copy of the class "exciton", + used here for test purposes only Linear chain of electronic sites with periodic boundaries Model consists of site energies (alpha) and nearest neighbor coupling (beta) From 96f2895247629505478ec54a93187fad45e9adeb Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Sat, 24 Jun 2023 11:10:17 +0200 Subject: [PATCH 06/25] Undo some of yesterday's changes because attributes "classical" and "bipartite" already existed --- wave_train/hamilton/coupled.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/wave_train/hamilton/coupled.py b/wave_train/hamilton/coupled.py index 3d2c581..3323735 100644 --- a/wave_train/hamilton/coupled.py +++ b/wave_train/hamilton/coupled.py @@ -76,8 +76,6 @@ def __init__(self, n_site, periodic, homogen, # Construct object of parent class Chain.__init__(self, n_site, periodic, homogen) - self.classical = True - self.bipartite = True # Construct objects of phonon and exciton class self.ex = Exciton(n_site, periodic, homogen, alpha, beta, eta) From 8abdb23ee6fa589dfc3aa22df1a8dac14e953097 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Sat, 24 Jun 2023 18:29:25 +0200 Subject: [PATCH 07/25] Redo previous changes of the import statements --- test_scripts/XYZ/tise_1.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test_scripts/XYZ/tise_1.py b/test_scripts/XYZ/tise_1.py index 30486c8..d827d65 100644 --- a/test_scripts/XYZ/tise_1.py +++ b/test_scripts/XYZ/tise_1.py @@ -1,6 +1,6 @@ -import wave_train.hamilton.xyz -import wave_train.dynamics.tise -import wave_train.io.logging +from wave_train.hamilton.xyz import XYZ +from wave_train.dynamics.tise import TISE +from wave_train.io.logging import TeeLogger from os.path import basename, splitext def xyz_tise(batch_mode): @@ -12,10 +12,10 @@ def xyz_tise(batch_mode): # class for logging to both console and logfile logger = None if not batch_mode: - logger = wave_train.io.logging.TeeLogger(log_file=my_file + ".log") + logger = TeeLogger(log_file=my_file + ".log") # Set up the excitonic Hamiltonian for a chain - hamilton = wave_train.hamilton.exciton.XYZ( + hamilton = XYZ( n_site=6, # number of sites periodic=True, # periodic boundary conditions homogen=True, # homogeneous chain/ring @@ -33,7 +33,7 @@ def xyz_tise(batch_mode): # Set up TISE solver # n_levels = n_site+1 # ground and singly excited # n_levels = n_site*(n_site+1)//2+1 # ground and singly and doubly excited - dynamics = wave_train.dynamics.tise.TISE( + dynamics = TISE( hamilton=hamilton, # choice of Hamiltonian, see above n_levels=16, # number of energy levels to be calculated solver='als', # choice of eigensolver for the full system From 8eed11ba234b10f08856fe77e58f1b5bd7a924a7 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 14 Jul 2023 11:49:00 +0200 Subject: [PATCH 08/25] Removing XYZ test files --- test_scripts/XYZ/tise_1.py | 75 -------------------------------------- 1 file changed, 75 deletions(-) delete mode 100644 test_scripts/XYZ/tise_1.py diff --git a/test_scripts/XYZ/tise_1.py b/test_scripts/XYZ/tise_1.py deleted file mode 100644 index d827d65..0000000 --- a/test_scripts/XYZ/tise_1.py +++ /dev/null @@ -1,75 +0,0 @@ -from wave_train.hamilton.xyz import XYZ -from wave_train.dynamics.tise import TISE -from wave_train.io.logging import TeeLogger -from os.path import basename, splitext - -def xyz_tise(batch_mode): - # Detect name of this script file (without extension) - base_name = basename(__file__) - my_file = splitext(base_name)[0] - - # logging instance: will be initialized with - # class for logging to both console and logfile - logger = None - if not batch_mode: - logger = TeeLogger(log_file=my_file + ".log") - - # Set up the excitonic Hamiltonian for a chain - hamilton = XYZ( - n_site=6, # number of sites - periodic=True, # periodic boundary conditions - homogen=True, # homogeneous chain/ring - alpha=0.1, # excitonic site energy - beta=-0.01, # coupling strength (NN) - eta=0.0 # constant energy offset - ) - - # Set up TT representation of the Hamiltonian - hamilton.get_TT( - n_basis=2, # size of electronic basis set - qtt=False # using quantized TT format - ) - - # Set up TISE solver - # n_levels = n_site+1 # ground and singly excited - # n_levels = n_site*(n_site+1)//2+1 # ground and singly and doubly excited - dynamics = TISE( - hamilton=hamilton, # choice of Hamiltonian, see above - n_levels=16, # number of energy levels to be calculated - solver='als', # choice of eigensolver for the full system - eigen='eig', # choice of eigensolver for the micro systems - ranks=15, # rank of initial guess for ALS - repeats=20, # number of sweeps in eigensolver scheme - conv_eps=1e-8, # threshold for detecting convergence of the eigenvalue - e_est=0.0, # estimation: eigenvalues closest to this number - e_min=0.05, # lower end of energy plot axis (if exact energies not available!) - e_max=0.15, # upper end of energy plot axis (if exact energies not available!) - save_file=my_file+'.pic', # if not None, generated data will be saved to this file - load_file=None, # if not None, reference data will be loaded from this file - compare=None # type of comparison with reference data - ) - - # Batch mode - if batch_mode: - dynamics.solve() # Solve TISE *without* visualization - - # Interactive mode: Setup animated visualization - else: - from wave_train.graphics.factory import VisualTISE - graphics = VisualTISE( - dynamics=dynamics, # choice of dynamics (EoM), see above - plot_type='QuantNumbers', # select your favorite plot type - plot_expect=True, # toggle plotting of expectation values - figure_pos=(100, 50), # specifying position (x,y) of upper left of figure [in pixels] - figure_size=(1050, 450), # specifying size (w,h) of figure [in pixels] - image_file=my_file+'.png', # if not None, image (last frame) will be written to this file - movie_file=my_file+'.mp4', # if not None, animation will be written to this file - snapshots=False, # each snapshot is saved - frame_rate=1, # frames per second in mp4 animation file - plot_style={} # additional plot style information - ).create() - graphics.solve() # Solve TISE *with* visualization - - -if __name__ == '__main__': - xyz_tise(batch_mode=False) From b6f80890bdf4313c0f8207b42c0beae525085243 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 14 Jul 2023 16:50:53 +0200 Subject: [PATCH 09/25] Removing xyz.py test file --- wave_train/hamilton/xyz.py | 308 ------------------------------------- 1 file changed, 308 deletions(-) delete mode 100644 wave_train/hamilton/xyz.py diff --git a/wave_train/hamilton/xyz.py b/wave_train/hamilton/xyz.py deleted file mode 100644 index 6dd8dec..0000000 --- a/wave_train/hamilton/xyz.py +++ /dev/null @@ -1,308 +0,0 @@ -import sys -import numpy as np -from wave_train.hamilton.chain import Chain - - -class XYZ(Chain): - - """ - Exciton dynamics in one dimension - --------------------------------- - - This is a copy of the class "exciton", - used here for test purposes only - Linear chain of electronic sites with periodic boundaries - Model consists of site energies (alpha) and nearest neighbor coupling (beta) - - Parent class is chain. - """ - - def __init__(self, n_site, periodic, homogen, alpha=1, beta=1, eta=0): - """ - Parameters: - ----------- - n_site: int - Length of the chain - periodic: boolean - Periodic boundary conditions - homogen: boolean - Homogeneous chain/ring - alpha: float - excitonic site energy - beta: float - coupling strength between nearest neighbours - eta: float - constant energy offset - """ - - # Initialize object of parent class - Chain.__init__(self, n_site, periodic, homogen) - self.classical = False - self.bipartite = False - - # Parameters of a homogeneous chain/ring are converted from scalar to vector - if homogen: - if isinstance(alpha,list): - sys.exit("Wrong input of alpha parameter: should be a scalar") - self.alpha = np.repeat(alpha, n_site) - - if isinstance(beta,list): - sys.exit("Wrong input of beta parameters: should be a scalar") - if periodic: - self.beta = np.repeat(beta, n_site) - else: - self.beta = np.repeat(beta, n_site-1) - - # Parameters of an inhomogeneous chain/ring should be given as vectors - else: - if len(alpha) != n_site: - sys.exit("Inconsistent length of vector of alpha parameters") - self.alpha = alpha - - if periodic: - if len(beta) != n_site: - sys.exit("Inconsistent length of vector of beta parameters") - else: - if len(beta) != n_site - 1: - sys.exit("Inconsistent length of vector of beta parameters") - self.beta = beta - - # Constant energy offset - self.eta = eta - - # Print output from superclass and from this class - print(super().__str__()) - print(self) - - # Useful for savemat|pickle output: Name and classification of this class - self.name = self.__class__.__name__ - self.bipartite = False # Single class of (quasi-) particles - self.classical = False # Classical approximation not supported - - def __str__(self): - if self.homogen: - alpha = self.alpha[0] - beta = self.beta[0] - else: - alpha = self.alpha - beta = self.beta - eta = self.eta - exciton_str = """ ------------------------- -Hamiltonian for EXCITONS ------------------------- - -Excitonic site energy (alpha) = {} -NN coupling strength (beta) = {} -Constant energy offset (eta) = {} - """.format(alpha, beta, eta) - - return exciton_str - - def get_2Q(self, n_dim=2): - """" - Second quantization: - Formulate Hamiltonian in terms of creation/annihiliation operators - """ - - # Call corresponding method from superclass Chain - super().get_2Q(n_dim) - - # Set up ground and first excited state - self.ground = np.zeros(n_dim) - self.ground[0] = 1 - self.excited = np.zeros(n_dim) - self.excited[1] = 1 - - def get_SLIM(self, n_dim=2): - """ - Set up SLIM matrices and operators of second quantization, - mainly for use in tensor train representation - - Parameters: - ----------- - n_dim: int - dimension of truncated(!) Hilbert space - - Returns - ------- - S: ndarray or list of ndarrays - single-site components of SLIM decomposition - L: ndarray or list of ndarrays - left two-site components of SLIM decomposition - I: ndarray or list of ndarrays - identity components of SLIM decomposition - M: ndarray or list of ndarrays - right two-site components of SLIM decomposition - - Reference: - Gelss/Klus/Matera/Schuette : Nearest-neighbor interaction systems in the tensor-train format - https://doi.org/10.1016/j.jcp.2017.04.007 - """ - - # Size of basis set - if n_dim != 2: - sys.exit(""" ---------------------------------- -Value of n_dim must be set to two ---------------------------------- - """) - self.n_dim = n_dim - - # Using second quantization - self.get_2Q (self.n_dim) - - # Set up S,L,I,M operators - S = self.qu_numbr - L_1 = self.raising - L_2 = self.lowering - I = self.identity - M_1 = self.lowering - M_2 = self.raising - - if self.homogen: - - # define ndarrays - self.S = self.alpha[0] * S + self.eta * I / self.n_site - self.L = self.beta[0] * np.stack((L_1, L_2), axis=-1) - self.I = I - self.M = np.stack((M_1, M_2)) - - else: - - # define lists of ndarrays - self.S = [None] * self.n_site - self.L = [None] * self.n_site - self.I = [None] * self.n_site - self.M = [None] * self.n_site - - # insert SLIM components - for i in range(self.n_site): # all sites (S) - self.S[i] = self.alpha[i] * S + self.eta * I / self.n_site - - for i in range(self.n_site - 1): # all sites but last (L) - self.L[i] = self.beta[i] * np.stack((L_1, L_2), axis=-1) - # self.L[i][:, :, 0] = self.beta[i] * L_1 - # self.L[i][:, :, 1] = self.beta[i] * L_2 - - for i in range(self.n_site): # all sites (I) - self.I[i] = I - - for i in range(1, self.n_site): # all sites but first (M) - self.M[i] = np.stack((M_1, M_2)) - # self.M[i][0, :, :] = M_1 - # self.M[i][1, :, :] = M_2 - - if self.periodic: # coupling last site (L) to first site (M) - self.L[self.n_site - 1] = self.beta[self.n_site - 1] * np.stack((L_1, L_2), axis=-1) - self.M[0] = np.stack((M_1, M_2)) - - print(""" -------------------------- -Setting up SLIM operators -------------------------- - -Size of excitonic basis set = {} - """.format(self.n_dim)) - - if self.homogen: - print('shape of S = ', self.S.shape) - print('shape of L = ', self.L.shape) - print('shape of I = ', self.I.shape) - print('shape of M = ', self.M.shape) - - def get_exact (self, n_levels): - """Compute analytic (exact) phononic energy levels for homogeneous systems, - i.e. analytic solutions of the underlying TISE. - So far, only up to two quanta of excitation - - Parameters: - ----------- - n_levels: int - number of exact energies to be computed - - Return: - ------ - energies: list of floats - list of exact energy levels - - Reference: - ---------- - Wikipedia article : Hueckel method - https://en.wikipedia.org/wiki/H%C3%BCckel_method - Z Hu, G. S. Engels, S. Kais : J. Chem. Phys. 148, 204307 (2018) - https://doi.org/10.1063/1.5026116 - """ - - # This method works only for homogeneous chains|rings - if not self.homogen: - return - - # Get necessary parameters - n_site = self.n_site - alpha = self.alpha[0] - beta = self.beta[0] - eta = self.eta - - # index j is an even integer running from 0 to 2N-2 - j = np.linspace(0, 2*n_site-2, n_site) - - # ground state: zero quanta of excitation - energies = eta * np.ones(1) - n_excite = 0 - - if self.periodic: # Cyclic systems (rings) - - # states with one quantum of excitation - if n_levels>1: - n_excite = 1 - ak = np.pi * j / n_site # wavenumber for lattice constant a=1 - ec = alpha + 2 * beta * np.cos(ak) # component energy - energies = np.append(energies, ec + eta) - - # states with two quanta of excitation - if n_levels>1+n_site: - n_excite = 2 - ak = np.pi * (j+1) / n_site # wavenumber for lattice constant a=1 - ec = alpha + 2 * beta * np.cos(ak) # component energy - ee,ff = np.meshgrid(ec,ec) # coordinate matrices from coordinate vectors - em = ee + ff # summing up two component energies - for i in range(1,n_site): # loop over all superdiagonals - energies = np.append(energies, np.diag(em,i)+ eta ) - - # states with more than two quanta not (yet?) implemented - if n_levels>1+sum(range(n_site+1)): - print ('Wrong choice of num_levels : ' + str(n_levels)) - sys.exit(1) - - else: # Linear systems (chains) - - # states with one quantum of excitation - if n_levels>1: - n_excite = 1 - ak = np.pi * (j/2+1) / (n_site+1) - ec = alpha + 2 * beta * np.cos(ak) - energies = np.append(energies, ec) - - # states with more than one quantum not (yet?) implemented - if n_levels>1+n_site: - print ('Wrong choice of num_levels : ' + str(n_levels)) - sys.exit(1) - - # sorting the energies in ascending order - energies = np.sort(energies) - - # truncate unnecessray energy levels - energies = energies[:n_levels] - - print(""" ------------------------------------------------- -Calculating excitonic energy levels analytically ------------------------------------------------- - -Number of energy levels desired : {} -Maximum number of excitations : {} - """.format(len(energies), n_excite)) - - return energies - From 58d6a6aee0273d8963000811dfccc8e480171e4d Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 14 Jul 2023 17:06:48 +0200 Subject: [PATCH 10/25] Abandon all occurrences of ex(citon) and ph(onon) throughout the 'dynamics' and 'graphics' classes --- wave_train/dynamics/qu_cl_mech.py | 8 ++--- wave_train/dynamics/quant_mech.py | 55 ++++++++++++++++--------------- wave_train/dynamics/tdse.py | 6 ++-- wave_train/dynamics/tise.py | 2 +- wave_train/graphics/services.py | 44 ++++++++++++------------- 5 files changed, 59 insertions(+), 56 deletions(-) diff --git a/wave_train/dynamics/qu_cl_mech.py b/wave_train/dynamics/qu_cl_mech.py index b16d67d..3bdeff1 100644 --- a/wave_train/dynamics/qu_cl_mech.py +++ b/wave_train/dynamics/qu_cl_mech.py @@ -32,7 +32,7 @@ def __init__(self, hamilton): # Quantum subsystem self.norm = np.zeros(self.num_steps + 1) # norm of state vector self.auto = np.zeros(self.num_steps + 1, dtype=complex) # autocorrelation - self.ex_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for excitons + self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for excitons # Classical subsystem self.position = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position @@ -91,16 +91,16 @@ def observe(self, i, iterations=0): for j in range(self.hamilton.n_site): self.position[i, j] = self.pos[j] self.momentum[i, j] = self.mom[j] - self.ex_numbr[i, j] = np.abs(self.psi[j]) ** 2 + self.q1_numbr[i, j] = np.abs(self.psi[j]) ** 2 - print(str("%4d" % j) + ' | ' + str("%10f" % self.ex_numbr[i, j]) + ' | ' + str("%10f" % self.position[i, j]) + ' | ' + str( + print(str("%4d" % j) + ' | ' + str("%10f" % self.q1_numbr[i, j]) + ' | ' + str("%10f" % self.position[i, j]) + ' | ' + str( "%10f" % self.momentum[i, j]) ) # Footer of table with site-specific information print(43 * '-') print(' sum' + - ' | ' + str("%10f" % np.sum(self.ex_numbr[i, :])) + + ' | ' + str("%10f" % np.sum(self.q1_numbr[i, :])) + ' | ' + str("%10f" % np.sum(self.position[i, :])) + ' | ' + str("%10f" % np.sum(self.momentum[i, :])) ) print (' ') diff --git a/wave_train/dynamics/quant_mech.py b/wave_train/dynamics/quant_mech.py index e09d7ae..10c2b81 100644 --- a/wave_train/dynamics/quant_mech.py +++ b/wave_train/dynamics/quant_mech.py @@ -46,8 +46,8 @@ def __init__(self, hamilton): self.unc_prod = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # uncertainty product if self.hamilton.bipartite: - self.ex_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for excitons - self.ph_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for phonons + self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 1st sub-system + self.q2_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 2nd sub-system else: self.qu_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number @@ -129,22 +129,23 @@ def observe(self, i, iterations=0): # Header of table with site-specific information print(' ') if self.hamilton.bipartite: - print('site | qu_n (ex) | qu_n (ph) | position | momentum | Dx (ph) | Dp (ph) | DxDp (ph)') + print('site | qu_n (1) | qu_n (2) | pos. (2) | mom. (2) | Dx (2) | Dp (2) | DxDp (2)') print(96 * '-') - elif self.hamilton.classical: - print('site | qu_number | position | momentum | Dx (ph) | Dp (ph) | DxDp (ph)') - print(82 * '-') else: - print( - 'site | qu_number') - print(17 * '-') + if self.hamilton.classical: + print('site | qu_number | position | momentum | Dx | Dp | DxDp') + print(82 * '-') + else: + print( + 'site | qu_number') + print(17 * '-') # Entries of table with site-specific information for j in range(self.hamilton.n_site): rho = self.rho_site[i,j] if self.hamilton.bipartite: - self.ex_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.ex_numbr), tol=TOLERANCE) - self.ph_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.ph_numbr), tol=TOLERANCE) + self.q1_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.q1_numbr), tol=TOLERANCE) + self.q2_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.q2_numbr), tol=TOLERANCE) else: self.qu_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.qu_numbr), tol=TOLERANCE) if self.hamilton.classical: @@ -157,29 +158,31 @@ def observe(self, i, iterations=0): self.unc_prod[i, j] = self.pos_wide[i,j] * self.mom_wide[i,j] # uncertainty product if self.hamilton.bipartite: - print(str("%4d" % j) + ' | ' + str("%10f" % self.ex_numbr[i, j]) + ' | ' + str( - "%10f" % self.ph_numbr[i, j]) + ' | ' + str("%10f" % self.position[i, j]) + ' | ' + str( + print(str("%4d" % j) + ' | ' + str("%10f" % self.q1_numbr[i, j]) + ' | ' + str( + "%10f" % self.q2_numbr[i, j]) + ' | ' + str("%10f" % self.position[i, j]) + ' | ' + str( "%10f" % self.momentum[i, j]) + ' | ' + str("%10f" % self.pos_wide[i, j]) + " | " + \ str("%10f" % self.mom_wide[i, j]) + " | " + str("%10f") % self.unc_prod[i, j] ) - elif self.hamilton.classical: - print(str("%4d" % j) + ' |' + str("%10f" % self.qu_numbr[i, j]) + ' | ' + str( - "%10f" % self.position[i, j]) + ' | ' + str("%10f" % self.momentum[i, j]) + ' | ' + \ - str("%10f" % self.pos_wide[i, j]) + ' | ' + str("%10f" % self.mom_wide[i, j]) + " | " + \ - str("%10f" % self.unc_prod[i, j])) else: - print(str("%4d" % j) + ' |' + str("%10f" % self.qu_numbr[i, j])) + if self.hamilton.classical: + print(str("%4d" % j) + ' |' + str("%10f" % self.qu_numbr[i, j]) + ' | ' + str( + "%10f" % self.position[i, j]) + ' | ' + str("%10f" % self.momentum[i, j]) + ' | ' + \ + str("%10f" % self.pos_wide[i, j]) + ' | ' + str("%10f" % self.mom_wide[i, j]) + " | " + \ + str("%10f" % self.unc_prod[i, j])) + else: + print(str("%4d" % j) + ' |' + str("%10f" % self.qu_numbr[i, j])) # Footer of table with site-specific information if self.hamilton.bipartite: print(96 * '-') - print(' sum | ' + str("%10f" % np.sum(self.ex_numbr[i, :])) + ' | ' + str( - "%10f" % np.sum(self.ph_numbr[i, :]))) - elif self.hamilton.classical: - print(82 * '-') - print(' sum |' + str("%10f" % np.sum(self.qu_numbr[i, :]))) + print(' sum | ' + str("%10f" % np.sum(self.q1_numbr[i, :])) + ' | ' + str( + "%10f" % np.sum(self.q2_numbr[i, :]))) else: - print(17 * '-') - print(' sum |' + str("%10f" % np.sum(self.qu_numbr[i, :]))) + if self.hamilton.classical: + print(82 * '-') + print(' sum |' + str("%10f" % np.sum(self.qu_numbr[i, :]))) + else: + print(17 * '-') + print(' sum |' + str("%10f" % np.sum(self.qu_numbr[i, :]))) print(' ') # RMSD of positions from reference solution (if available) diff --git a/wave_train/dynamics/tdse.py b/wave_train/dynamics/tdse.py index 062b6ae..3ac1105 100644 --- a/wave_train/dynamics/tdse.py +++ b/wave_train/dynamics/tdse.py @@ -23,7 +23,7 @@ def __init__(self, hamilton, num_steps, step_size, sub_steps=1, save_file=None, load_file=None, compare=None): """ hamilton: instance of physical object (quantum Hamiltonian) - Either one Exciton, Phonon, Coupled or other classes + Either one of Exciton, Phonon, Coupled or other classes num_steps: int number of (main) time steps step_size: float @@ -288,8 +288,8 @@ def start_solve(self): self.e_min = self.e_init*0.9 self.e_max = self.e_init*1.1 - # Analytic solutions in terms of Bessel functions for excitons on an homogenous chain - if self.hamilton.name == 'Exciton' and self.hamilton.homogen and not self.hamilton.periodic: + # Analytic solutions in terms of Bessel functions: for linear homogenous 2-state systems only! + if self.hamilton.n_dim == 2 and self.hamilton.homogen and not self.hamilton.periodic: # If only a single site is initially excited if hasattr(self, 'localized'): diff --git a/wave_train/dynamics/tise.py b/wave_train/dynamics/tise.py index c58b0d3..ffbe533 100644 --- a/wave_train/dynamics/tise.py +++ b/wave_train/dynamics/tise.py @@ -22,7 +22,7 @@ def __init__(self, hamilton, n_levels, Parameters: ----------- hamilton: instance of physical object (quantum Hamiltonian) - Either one Exciton, Phonon, Coupled or other classes + Either one of Exciton, Phonon, Coupled or other classes n_levels: int number of solutions (energy levels) of the TISE solver: string, optional diff --git a/wave_train/graphics/services.py b/wave_train/graphics/services.py index c81c476..3293b1d 100644 --- a/wave_train/graphics/services.py +++ b/wave_train/graphics/services.py @@ -314,15 +314,15 @@ def configure_quant_numbers2_basic(figure, dynamics, outer_grid=None, style=figu axes_number.set_yticks(np.linspace(0, 2, 5)) axes_number.set_ylabel("Quantum numbers") - # generate initial bars of height zero for excitonic and phononic system - bars_exc = np.arange(dynamics.hamilton.n_site) - 0.1 - bars_pho = np.arange(dynamics.hamilton.n_site) + 0.1 + # generate initial bars of height zero for the two sub-systems + bars_s1 = np.arange(dynamics.hamilton.n_site) - 0.1 + bars_s2 = np.arange(dynamics.hamilton.n_site) + 0.1 heights = np.zeros(dynamics.hamilton.n_site) - # first n_site rectangles belong to excitonic quantum numbers only, - # second half to phononic quantum numbers only - axes_number.bar(bars_exc, heights, width=0.2, color='green', label="Quantum") - axes_number.bar(bars_pho, heights, width=0.2, color='orange', label="Classical") + # first n_site rectangles belong to 1st sub-system quantum numbers only, + # second half to 2nd sub-system quantum numbers only + axes_number.bar(bars_s1, heights, width=0.2, color='green', label="Quantum") + axes_number.bar(bars_s2, heights, width=0.2, color='orange', label="Classical") axes_number.legend() apply_styles([axes_number], style) @@ -657,7 +657,7 @@ def configure_quant_displace2_expect_qcmd(figure, dynamics, style=figure_style): Quantum numbers for each site are retrieved from the (reduced) density matrices and the number operators using the trace formula for expectation values. - These operations are done separately for excitonic and phononic densities + These operations are done separately for the two sub-sytems' densities """ outer_grid = gridspec.GridSpec(3, 2, figure=figure) configure_quant_displace2_basic(figure, dynamics, outer_grid, style=style) @@ -745,13 +745,13 @@ def update_quant_numbers2_basic(i, figure, dynamics, writer, saving): dynamics.update_solve(i) # get all rectangle instances of the first plot - # first half of array corresponds to excitonic values - # second half of array corresponds to phononic values + # first half of array corresponds to 1st sub-system values + # second half of array corresponds to 2nd sub-system values rectangles = figure.axes[0].patches for k in range(dynamics.hamilton.n_site): - rectangles[k].set_height(dynamics.ex_numbr[i, k]) - rectangles[k + dynamics.hamilton.n_site].set_height(dynamics.ph_numbr[i, k]) + rectangles[k].set_height(dynamics.q1_numbr[i, k]) + rectangles[k + dynamics.hamilton.n_site].set_height(dynamics.q2_numbr[i, k]) figure.suptitle(dynamics.head[i]) @@ -765,7 +765,7 @@ def update_quant_displace2_basic(i, figure, dynamics, writer, saving): # split densities and get expectation values of number operators for j in range(dynamics.hamilton.n_site): - rectangles[j].set_height(dynamics.ex_numbr[i, j]) + rectangles[j].set_height(dynamics.q1_numbr[i, j]) height = scale * dynamics.position[i, j] rectangles[j + dynamics.hamilton.n_site].set_height(height) @@ -777,17 +777,17 @@ def update_positions2_basic(i, figure, dynamics, writer, saving): dynamics.update_solve(i) pos = scaling() * dynamics.position[i] - nex = dynamics.ex_numbr[i] + n_1 = dynamics.q1_numbr[i] - # all even indices correspond to phononic system - phononic = np.asarray(figure.axes[0].lines)[0::2] - # all odd indices correspond to excitonic system - excitonic = np.asarray(figure.axes[0].lines)[1::2] + # all even indices correspond to 2nd sub-system + s2 = np.asarray(figure.axes[0].lines)[0::2] + # all odd indices correspond to 1st sub-system + s1 = np.asarray(figure.axes[0].lines)[1::2] - for j, (line_ph, line_ex) in enumerate(zip(phononic, excitonic)): - line_ph.set_data([j, j + pos[j]], [0, 0]) - # move excitonic quantum number with point of displacement - line_ex.set_data([j + pos[j], j + pos[j]], [0, nex[j]]) + for j, (line_2, line_1) in enumerate(zip(s2, s1)): + line_2.set_data([j, j + pos[j]], [0, 0]) + # move 1st sub-system quantum number with point of displacement + line_1.set_data([j + pos[j], j + pos[j]], [0, n_1[j]]) figure.suptitle(dynamics.head[i]) From 65b1c8fc031b7a5299ed9e52d796f86fd36fb0f7 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 14 Jul 2023 17:47:01 +0200 Subject: [PATCH 11/25] really tiny change --- wave_train/dynamics/qu_cl_mech.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wave_train/dynamics/qu_cl_mech.py b/wave_train/dynamics/qu_cl_mech.py index 3bdeff1..427dbb2 100644 --- a/wave_train/dynamics/qu_cl_mech.py +++ b/wave_train/dynamics/qu_cl_mech.py @@ -32,7 +32,7 @@ def __init__(self, hamilton): # Quantum subsystem self.norm = np.zeros(self.num_steps + 1) # norm of state vector self.auto = np.zeros(self.num_steps + 1, dtype=complex) # autocorrelation - self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for excitons + self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 1st sub-system # Classical subsystem self.position = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position From 19e9cb285ca02cd4921a5e25e863c7e3a4df3e64 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 14 Jul 2023 18:02:14 +0200 Subject: [PATCH 12/25] really tiny change --- wave_train/graphics/animation.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/wave_train/graphics/animation.py b/wave_train/graphics/animation.py index 9400e0a..4a874b1 100644 --- a/wave_train/graphics/animation.py +++ b/wave_train/graphics/animation.py @@ -14,7 +14,7 @@ class Animation: solve the underlying problem of the FuncAnimation class, namely that creation of animation output and visualization are executed in order and not in parallel. - 1) If image_file is not None, last snaphot is saved, but previous snapshots are removed. + 1) If image_file is not None, last snapshot is saved, but previous snapshots are removed. 2) If image_file is not None and snapshot is True, all snapshots are saved. If image_file == None and snapshot == True, a default snapshot name will be assigned 3) If image_file is not None, movie_file is not None, and snapshot is True, @@ -34,7 +34,7 @@ class Animation: "frames", # frames in the animation "pause", # pause between frames "frame_rate", # number of frames per second - "frame_count" # counter of frames/snaphots in animation + "frame_count" # counter of frames/snapshots in animation ] movie_formats = ['mp4'] @@ -87,10 +87,10 @@ def __init__(self, dynamics, movie_file, image_file, snapshots, figure, update, if self.image_file is None: if snapshots and self.movie_file is not None: - self.image_file = "snaphot.png" + self.image_file = "snapshot.png" logging.warning("Name for snapshots was not provided. Defaulting to 'snapshot.png'") elif self.movie_file is not None: - # indicate that all snaphots must be removed + # indicate that all snapshots must be removed self.image_file = "snapshot_removable.png" else: self.image_file = "" From e5b92a53b1f98c1c99307af858ba8b4e01455fda Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 14 Jul 2023 19:26:52 +0200 Subject: [PATCH 13/25] Calculate classical variables only where needed --- wave_train/dynamics/quant_mech.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/wave_train/dynamics/quant_mech.py b/wave_train/dynamics/quant_mech.py index 10c2b81..6087683 100644 --- a/wave_train/dynamics/quant_mech.py +++ b/wave_train/dynamics/quant_mech.py @@ -37,13 +37,6 @@ def __init__(self, hamilton): # More observables (from reduced density matrices) self.rho_site = np.zeros((self.num_steps + 1, self.hamilton.n_site, self.hamilton.n_dim, self.hamilton.n_dim), dtype=complex) - self.position = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position - self.momentum = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum - self.pos_squa = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position squared - self.mom_squa = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum uncertainty - self.pos_wide = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position uncertainty - self.mom_wide = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum uncertainty - self.unc_prod = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # uncertainty product if self.hamilton.bipartite: self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 1st sub-system @@ -51,6 +44,15 @@ def __init__(self, hamilton): else: self.qu_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number + if self.hamilton.classical: + self.position = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position + self.momentum = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum + self.pos_squa = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position squared + self.mom_squa = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum uncertainty + self.pos_wide = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position uncertainty + self.mom_wide = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum uncertainty + self.unc_prod = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # uncertainty product + def observe(self, i, iterations=0): """ Update and print expectation values of observables @@ -140,7 +142,7 @@ def observe(self, i, iterations=0): 'site | qu_number') print(17 * '-') - # Entries of table with site-specific information + # Entries of table with site-specific information for j in range(self.hamilton.n_site): rho = self.rho_site[i,j] if self.hamilton.bipartite: From edee6945821e74d9d206980edc818fb8a5c4897a Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Tue, 1 Aug 2023 15:38:16 +0200 Subject: [PATCH 14/25] Replaced 'se' by 's2' (for symmetric Euler aka SOD) throughout. But why does 'se' not throw an error??? --- test_scripts/Coupled/tdse_1.py | 2 +- test_scripts/Coupled/tdse_qe_2.py | 2 +- test_scripts/Exciton/tdse_1.py | 2 +- test_scripts/Exciton/tdse_qe_2.py | 2 +- test_scripts/Phonon/tdse_1.py | 2 +- test_scripts/Phonon/tdse_qe_2.py | 2 +- 6 files changed, 6 insertions(+), 6 deletions(-) diff --git a/test_scripts/Coupled/tdse_1.py b/test_scripts/Coupled/tdse_1.py index ec6e91c..9a43050 100644 --- a/test_scripts/Coupled/tdse_1.py +++ b/test_scripts/Coupled/tdse_1.py @@ -43,7 +43,7 @@ def coupled_tdse(batch_mode): num_steps=25, # number of main time steps step_size=0.01/hamilton.sig[0], # size of main time steps sub_steps=10, # number of sub steps - solver='sm', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + solver='sm', # can be 's2' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=12, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) diff --git a/test_scripts/Coupled/tdse_qe_2.py b/test_scripts/Coupled/tdse_qe_2.py index ec55e87..32ce4bb 100644 --- a/test_scripts/Coupled/tdse_qe_2.py +++ b/test_scripts/Coupled/tdse_qe_2.py @@ -43,7 +43,7 @@ def coupled_qe_2(batch_mode): num_steps=50, # number of main time steps step_size=20, # size of main time steps sub_steps=20, # number of sub steps - solver='sm', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + solver='sm', # can be 's2' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=12, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) diff --git a/test_scripts/Exciton/tdse_1.py b/test_scripts/Exciton/tdse_1.py index ef0b471..5bf0512 100644 --- a/test_scripts/Exciton/tdse_1.py +++ b/test_scripts/Exciton/tdse_1.py @@ -36,7 +36,7 @@ def exciton_tdse(batch_mode): num_steps=50, # number of main time steps step_size=20, # size of main time steps sub_steps=5, # number of sub steps - solver='sm', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + solver='sm', # can be 's2' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=8, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) diff --git a/test_scripts/Exciton/tdse_qe_2.py b/test_scripts/Exciton/tdse_qe_2.py index 49c73a6..e0cf4c0 100644 --- a/test_scripts/Exciton/tdse_qe_2.py +++ b/test_scripts/Exciton/tdse_qe_2.py @@ -36,7 +36,7 @@ def exciton_qe_2(batch_mode): num_steps=10, # number of main time steps step_size=20, # size of main time steps sub_steps=10, # number of sub steps - solver='sm', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + solver='sm', # can be 's2' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=8, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) diff --git a/test_scripts/Phonon/tdse_1.py b/test_scripts/Phonon/tdse_1.py index b5e9974..dabe7a5 100644 --- a/test_scripts/Phonon/tdse_1.py +++ b/test_scripts/Phonon/tdse_1.py @@ -36,7 +36,7 @@ def phonon_tdse(batch_mode): num_steps=50, # number of main time steps step_size=200, # size of main time steps sub_steps=50, # number of sub steps - solver='se', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + solver='s2', # can be 's2' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=15, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) diff --git a/test_scripts/Phonon/tdse_qe_2.py b/test_scripts/Phonon/tdse_qe_2.py index deb4982..15345e4 100644 --- a/test_scripts/Phonon/tdse_qe_2.py +++ b/test_scripts/Phonon/tdse_qe_2.py @@ -36,7 +36,7 @@ def phonon_qe_2(batch_mode): num_steps=50, # number of main time steps step_size=200, # size of main time steps sub_steps=50, # number of sub steps - solver='se', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + solver='s2', # can be 's2' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=15, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) From f68e8b6c633cbce53b208ef91863d36c912a7916 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Tue, 1 Aug 2023 16:58:36 +0200 Subject: [PATCH 15/25] Renamed ex_, ph_ into q1_, q2_. Now running OK! --- wave_train/hamilton/coupled.py | 50 +++++++++++++++++----------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/wave_train/hamilton/coupled.py b/wave_train/hamilton/coupled.py index 3323735..f64db79 100644 --- a/wave_train/hamilton/coupled.py +++ b/wave_train/hamilton/coupled.py @@ -207,16 +207,16 @@ def get_2Q(self, n_dims): self.mom_conv = self.ph.mom_conv # Excitonic Hamiltonian in terms of creation/annihilation operators - self.ex_raise = np.kron(self.ex.raising, self.ph.identity) # raising operator: excitons - self.ex_lower = np.kron(self.ex.lowering, self.ph.identity) # lowering operator: excitons - self.ex_numbr = np.kron(self.ex.qu_numbr, self.ph.identity) # number operator: excitons + self.q1_raise = np.kron(self.ex.raising, self.ph.identity) # raising operator: excitons + self.q1_lower = np.kron(self.ex.lowering, self.ph.identity) # lowering operator: excitons + self.q1_numbr = np.kron(self.ex.qu_numbr, self.ph.identity) # number operator: excitons # Phononic Hamiltonian in terms of creation/annihilation operators - self.ph_raise = np.kron(self.ex.identity, self.ph.raising) # raising operator: phonons - self.ph_lower = np.kron(self.ex.identity, self.ph.lowering) # lowering operator: phonons - self.ph_numbr = np.kron(self.ex.identity, self.ph.qu_numbr) # number operator: phonons - self.position = self.ph_raise + self.ph_lower # position operator: phonons - self.momentum = self.ph_raise - self.ph_lower # momentum operator: phonons + self.q2_raise = np.kron(self.ex.identity, self.ph.raising) # raising operator: phonons + self.q2_lower = np.kron(self.ex.identity, self.ph.lowering) # lowering operator: phonons + self.q2_numbr = np.kron(self.ex.identity, self.ph.qu_numbr) # number operator: phonons + self.position = self.q2_raise + self.q2_lower # position operator: phonons + self.momentum = self.q2_raise - self.q2_lower # momentum operator: phonons self.pos_squa = self.position @ self.position self.mom_squa = self.momentum @ self.momentum @@ -318,33 +318,33 @@ def get_SLIM(self, n_dims): sig_2[i] *= self.pos_conv[i] # Set up S,L,I,M operators - S_1 = self.ex_numbr # on site: excitons - S_2 = self.ph_numbr + self.identity / 2 # on site: phonons - S_3 = self.ex_numbr @ self.position # on site: exciton-phonon-tuning (matrix product!) + S_1 = self.q1_numbr # on site: excitons + S_2 = self.q2_numbr + self.identity / 2 # on site: phonons + S_3 = self.q1_numbr @ self.position # on site: exciton-phonon-tuning (matrix product!) - L_1 = self.ex_raise # NN interaction: exciton transfer - L_2 = self.ex_lower # NN interaction: exciton transfer + L_1 = self.q1_raise # NN interaction: exciton transfer + L_2 = self.q1_lower # NN interaction: exciton transfer L_3 = self.position # NN interaction: phonon coupling - L_4 = self.ex_numbr # NN interaction: ex-ph-tuning + L_4 = self.q1_numbr # NN interaction: ex-ph-tuning L_5 = self.position # NN interaction: ex-ph-tuning if tau_1[0] != 0: # NN interaction: ex-ph-coupling - L_6 = self.ex_raise - L_7 = self.ex_raise @ self.position - L_8 = self.ex_lower - L_9 = self.ex_lower @ self.position + L_6 = self.q1_raise + L_7 = self.q1_raise @ self.position + L_8 = self.q1_lower + L_9 = self.q1_lower @ self.position I = self.identity # identity operator - M_1 = self.ex_lower # NN interaction: exciton transfer - M_2 = self.ex_raise # NN interaction: exciton transfer + M_1 = self.q1_lower # NN interaction: exciton transfer + M_2 = self.q1_raise # NN interaction: exciton transfer M_3 = self.position # NN interaction: phonon coupling M_4 = self.position # NN interaction: ex-ph-tuning - M_5 = self.ex_numbr # NN interaction: ex-ph-tuning + M_5 = self.q1_numbr # NN interaction: ex-ph-tuning if tau_1[0] != 0: # NN interaction: ex-ph-coupling - M_6 = self.ex_lower @ self.position - M_7 = self.ex_lower - M_8 = self.ex_raise @ self.position - M_9 = self.ex_raise + M_6 = self.q1_lower @ self.position + M_7 = self.q1_lower + M_8 = self.q1_raise @ self.position + M_9 = self.q1_raise # Due to the use of the efficient frequencies nu_E, omg_E, # non-periodic phonon chains appear here as non-homogeneous From 52c1962b2853159e07c9f547c7cd86895b8aff64 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Wed, 2 Aug 2023 17:23:41 +0200 Subject: [PATCH 16/25] Implemented the three Pauli matrices, along with calculating/printing expectation values --- wave_train/dynamics/quant_mech.py | 44 ++++++++++++++++++------------- wave_train/dynamics/tdse.py | 2 -- wave_train/hamilton/chain.py | 13 +++++++-- 3 files changed, 37 insertions(+), 22 deletions(-) diff --git a/wave_train/dynamics/quant_mech.py b/wave_train/dynamics/quant_mech.py index 6087683..6810ce5 100644 --- a/wave_train/dynamics/quant_mech.py +++ b/wave_train/dynamics/quant_mech.py @@ -39,11 +39,13 @@ def __init__(self, hamilton): self.hamilton.n_dim, self.hamilton.n_dim), dtype=complex) if self.hamilton.bipartite: - self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 1st sub-system - self.q2_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 2nd sub-system + self.q1_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 1st sub-system + self.q2_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number for 2nd sub-system else: - self.qu_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number - + self.qu_numbr = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # quantum number + self.qu_sig_1 = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # Pauli matrix + self.qu_sig_2 = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # Pauli matrix + self.qu_sig_3 = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # Pauli matrix if self.hamilton.classical: self.position = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # position self.momentum = np.zeros((self.num_steps + 1, self.hamilton.n_site)) # momentum @@ -139,25 +141,28 @@ def observe(self, i, iterations=0): print(82 * '-') else: print( - 'site | qu_number') - print(17 * '-') + 'site | qu_number | sigma_1 | sigma_2 | sigma_3') + print(52 * '-') # Entries of table with site-specific information for j in range(self.hamilton.n_site): rho = self.rho_site[i,j] if self.hamilton.bipartite: - self.q1_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.q1_numbr), tol=TOLERANCE) - self.q2_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.q2_numbr), tol=TOLERANCE) + self.q1_numbr[i, j] = np.real_if_close(np.trace(rho @ self.hamilton.q1_numbr), tol=TOLERANCE) + self.q2_numbr[i, j] = np.real_if_close(np.trace(rho @ self.hamilton.q2_numbr), tol=TOLERANCE) else: - self.qu_numbr[i, j] = np.real_if_close( np.trace(rho @ self.hamilton.qu_numbr), tol=TOLERANCE) + self.qu_numbr[i, j] = np.real_if_close(np.trace(rho @ self.hamilton.qu_numbr), tol=TOLERANCE) + self.qu_sig_1[i, j] = np.real_if_close(np.trace(rho @ self.hamilton.qu_sig_1), tol=TOLERANCE) + self.qu_sig_2[i, j] = np.real_if_close(np.trace(rho @ self.hamilton.qu_sig_2), tol=TOLERANCE) + self.qu_sig_3[i, j] = np.real_if_close(np.trace(rho @ self.hamilton.qu_sig_3), tol=TOLERANCE) if self.hamilton.classical: - self.position[i, j] = np.real_if_close(self.hamilton.pos_conv[j] * np.trace(rho @ self.hamilton.position), tol=TOLERANCE) - self.momentum[i, j] = np.real_if_close(self.hamilton.mom_conv[j] * np.trace(rho @ self.hamilton.momentum), tol=TOLERANCE) - self.pos_squa[i, j] = np.real_if_close(self.hamilton.pos_conv[j]**2 * np.trace(rho @ self.hamilton.pos_squa), tol=TOLERANCE) - self.mom_squa[i, j] = np.real_if_close(self.hamilton.mom_conv[j]**2 * np.trace(rho @ self.hamilton.mom_squa), tol=TOLERANCE) - self.pos_wide[i, j] = np.real_if_close( np.sqrt(self.pos_squa[i,j] - self.position[i,j]**2), tol=TOLERANCE) - self.mom_wide[i, j] = np.real_if_close( np.sqrt(self.mom_squa[i,j] - self.momentum[i,j]**2), tol=TOLERANCE) - self.unc_prod[i, j] = self.pos_wide[i,j] * self.mom_wide[i,j] # uncertainty product + self.position[i, j] = np.real_if_close(self.hamilton.pos_conv[j] * np.trace(rho @ self.hamilton.position), tol=TOLERANCE) + self.momentum[i, j] = np.real_if_close(self.hamilton.mom_conv[j] * np.trace(rho @ self.hamilton.momentum), tol=TOLERANCE) + self.pos_squa[i, j] = np.real_if_close(self.hamilton.pos_conv[j]**2 * np.trace(rho @ self.hamilton.pos_squa), tol=TOLERANCE) + self.mom_squa[i, j] = np.real_if_close(self.hamilton.mom_conv[j]**2 * np.trace(rho @ self.hamilton.mom_squa), tol=TOLERANCE) + self.pos_wide[i, j] = np.real_if_close(np.sqrt(self.pos_squa[i,j] - self.position[i,j]**2), tol=TOLERANCE) + self.mom_wide[i, j] = np.real_if_close(np.sqrt(self.mom_squa[i,j] - self.momentum[i,j]**2), tol=TOLERANCE) + self.unc_prod[i, j] = self.pos_wide[i,j] * self.mom_wide[i,j] # uncertainty product if self.hamilton.bipartite: print(str("%4d" % j) + ' | ' + str("%10f" % self.q1_numbr[i, j]) + ' | ' + str( @@ -171,7 +176,10 @@ def observe(self, i, iterations=0): str("%10f" % self.pos_wide[i, j]) + ' | ' + str("%10f" % self.mom_wide[i, j]) + " | " + \ str("%10f" % self.unc_prod[i, j])) else: - print(str("%4d" % j) + ' |' + str("%10f" % self.qu_numbr[i, j])) + print(str("%4d" % j) + ' |' + str("%10f" % self.qu_numbr[i, j]) + \ + ' |' + str("%10f" % self.qu_sig_1[i, j]) + \ + ' |' + str("%10f" % self.qu_sig_2[i, j]) + \ + ' |' + str("%10f" % self.qu_sig_3[i, j])) # Footer of table with site-specific information if self.hamilton.bipartite: @@ -183,7 +191,7 @@ def observe(self, i, iterations=0): print(82 * '-') print(' sum |' + str("%10f" % np.sum(self.qu_numbr[i, :]))) else: - print(17 * '-') + print(52 * '-') print(' sum |' + str("%10f" % np.sum(self.qu_numbr[i, :]))) print(' ') diff --git a/wave_train/dynamics/tdse.py b/wave_train/dynamics/tdse.py index 3ac1105..658617b 100644 --- a/wave_train/dynamics/tdse.py +++ b/wave_train/dynamics/tdse.py @@ -433,11 +433,9 @@ def update_solve(self, i): initial_value=self.psi, step_size=self.sub_size, number_of_steps=self.sub_steps, threshold=self.threshold, max_rank=self.max_rank, normalize=self.normalize) - elif self.solver == 'vp': # time-dependent variational principle psi = ode.tdvp(1j*self.operator, initial_value=self.psi, step_size=self.sub_size, number_of_steps=self.sub_steps, normalize=self.normalize) - elif self.solver == 'qe': # Quasi-exact propagation self.psi = self.propagator @ self.psi diff --git a/wave_train/hamilton/chain.py b/wave_train/hamilton/chain.py index 56c8813..8041192 100644 --- a/wave_train/hamilton/chain.py +++ b/wave_train/hamilton/chain.py @@ -72,8 +72,17 @@ def get_2Q(self, n_dim): self.qu_numbr = self.raising @ self.lowering # number operator self.position = self.raising + self.lowering # position operator self.momentum = self.raising - self.lowering # momentum operator - self.pos_squa = self.position @ self.position - self.mom_squa = self.momentum @ self.momentum + self.pos_squa = self.position @ self.position # position squared + self.mom_squa = self.momentum @ self.momentum # momentum squared + self.qu_sig_1 = np.zeros((n_dim,n_dim)) # 1st Pauli matrix + self.qu_sig_1[0, 1] = 1 + self.qu_sig_1[1, 0] = 1 + self.qu_sig_2 = np.zeros((n_dim,n_dim), dtype=complex) # 2nd Pauli matrix + self.qu_sig_2[0, 1] = -1j + self.qu_sig_2[1, 0] = 1j + self.qu_sig_3 = np.zeros((n_dim,n_dim)) # 3rd Pauli matrix + self.qu_sig_3[0, 0] = 1 + self.qu_sig_3[1, 1] = -1 def get_TT(self, n_basis, qtt=False): """ From ecf62b3c78a14661d4591fc28a19372ac6b0d6d6 Mon Sep 17 00:00:00 2001 From: PGelss Date: Fri, 1 Sep 2023 14:55:20 +0200 Subject: [PATCH 17/25] added boson.py for mapping of bath degrees of freedom to a 1d chain --- wave_train/hamilton/boson.py | 162 +++++++++++++++++++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 wave_train/hamilton/boson.py diff --git a/wave_train/hamilton/boson.py b/wave_train/hamilton/boson.py new file mode 100644 index 0000000..7e9b359 --- /dev/null +++ b/wave_train/hamilton/boson.py @@ -0,0 +1,162 @@ +import sys +import numpy as np +from wave_train.hamilton.chain import Chain + + +class Boson(Chain): + + """ + Boson ... + --------------------------------- + + ... + + Parent class is chain. + """ + + def __init__(self, n_site, eta, s, omega_c, omega_0): + """ + Parameters: + ----------- + n_site: int + Length of the chain + eta: float + coupling between TLS and first bath site + s: float + type of spectral density function: s<1 sub-ohmic, s=1 ohmic, s>1 super-ohmic + omega_c: float + cut-off frequency of spectral density function + omega_0: float + eigenfrequency of the TLS + """ + + # Initialize object of parent class + Chain.__init__(self, n_site, False, False) + + self.eta = eta + self.s = s + self.omega_c = omega_c + self.omega_0 = omega_0 + + # Print output from superclass and from this class + print(super().__str__()) + print(self) + + # Useful for savemat|pickle output: Name and classification of this class + self.name = self.__class__.__name__ + self.bipartite = False # Single class of (quasi-) particles + self.classical = False # Classical approximation not supported + + def __str__(self): + eta = self.eta + s = self.s + omega_c = self.omega_c + omega_0 = self.omega_0 + + boson_str = """ +---------------------- +Hamiltonian for BOSONS +---------------------- + +Coupling between TLS and first bath site (eta) = {} +Type of spectral density function (s) = {} +Cut-off frequency of spectral density function (omega_c) = {} +Eigenfrequency of the TLS (omega_0) = {} + """.format(eta, s, omega_c, omega_0) + + return boson_str + + def get_2Q(self, n_dim=2): + """" + Second quantization: + Formulate Hamiltonian in terms of creation/annihiliation operators + """ + + # Call corresponding method from superclass Chain + super().get_2Q(n_dim) + + # Set up ground and first excited state + self.ground = np.zeros(n_dim) + self.ground[0] = 1 + self.excited = np.zeros(n_dim) + self.excited[1] = 1 + self.cat = np.zeros(n_dim) + self.cat[0] = np.sqrt(0.5) + self.cat[1] = np.sqrt(0.5) + + + def get_SLIM(self, n_dim=2): + """ + Set up SLIM matrices and operators of second quantization, + mainly for use in tensor train representation + + Parameters: + ----------- + n_dim: int + dimension of truncated(!) Hilbert space + + Returns + ------- + S: ndarray or list of ndarrays + single-site components of SLIM decomposition + L: ndarray or list of ndarrays + left two-site components of SLIM decomposition + I: ndarray or list of ndarrays + identity components of SLIM decomposition + M: ndarray or list of ndarrays + right two-site components of SLIM decomposition + + Reference: + Gelss/Klus/Matera/Schuette : Nearest-neighbor interaction systems in the tensor-train format + https://doi.org/10.1016/j.jcp.2017.04.007 + """ + + self.n_dim = 3 + + # Using second quantization + self.get_2Q (self.n_dim) + + sigmax = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]) + sigmay = np.array([[0, -1j, 0], [1j, 0, 0], [0, 0, 0]]) + sigmaz = np.array([[1, 0, 0], [0, -1, 0], [0, 0, 0]]) + #nbtotal = self.qu_numbr# np.diag([0, 1, 2]) + #b = self.lowering#np.diag([1,np.sqrt(2)], k=1) + #bdagger = self.raising# np.diag([1,np.sqrt(2)], k=-1) + #Q = self.position#b + bdagger + + # bond lists + tls_chain_bond = [self.eta] + [0]*(self.n_site-2) + chain_bond = [0] + [self.omega_c*(1+i)*(1+i+self.s)/((2+2*i+self.s)*(3+2*i+self.s))*np.sqrt((3+2*i+self.s)/(1+2*i+self.s)) for i in range(self.n_site-2)] + tls_onsite = [self.omega_0] + [0]*(self.n_site-1) + cavity_onsite = [0] + [self.omega_c/2*(1 + self.s**2/((self.s+2*i)*(2+self.s+2*i))) for i in range(self.n_site-1)] + + + # define lists of ndarrays + self.S = [self.omega_0*sigmaz] + [cavity_onsite[i]*self.qu_numbr for i in range(1,self.n_site)] + self.L = [np.stack((self.eta*sigmax, self.eta*self.position), axis=-1)] + [np.stack((chain_bond[i]*self.lowering,chain_bond[i]*self.raising), axis=-1) for i in range(1,self.n_site-1)] + [None] + self.I = [np.eye(3) for i in range(self.n_site)] + self.M = [None, np.stack((self.position,sigmax))] + [np.stack((self.raising,self.lowering)) for i in range(2,self.n_site)] + + + +# OBSERVABLES +# +# given a quantum state psi as TT, do the following: +# +# sigmax = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0]]) +# sigmay = np.array([[0, -1j, 0], [1j, 0, 0], [0, 0, 0]]) +# sigmaz = np.array([[1, 0, 0], [0, -1, 0], [0, 0, 0]]) +# +# exp_vals = np.zeros([3, psi.order]) +# for i in range(psi.order): +# psi_tmp = psi.copy() +# psi_tmp.cores[i] = np.einsum('ij,kjlm->kilm', sigmax, psi_tmp.cores[i]) +# exp_vals[0, i] = psi.transpose(conjugate=True)@psi_tmp +# psi_tmp = psi.copy() +# psi_tmp.cores[i] = np.einsum('ij,kjlm->kilm', sigmay, psi_tmp.cores[i]) +# exp_vals[1, i] = psi.transpose(conjugate=True)@psi_tmp +# psi_tmp = psi.copy() +# psi_tmp.cores[i] = np.einsum('ij,kjlm->kilm', sigmaz, psi_tmp.cores[i]) +# exp_vals[2, i] = psi.transpose(conjugate=True)@psi_tmp + + From b2251c34f119f06892fd8d5403f1db696b829dd0 Mon Sep 17 00:00:00 2001 From: PGelss Date: Fri, 22 Sep 2023 11:09:25 +0200 Subject: [PATCH 18/25] added Krylov method --- wave_train/dynamics/tdse.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/wave_train/dynamics/tdse.py b/wave_train/dynamics/tdse.py index 658617b..26234c9 100644 --- a/wave_train/dynamics/tdse.py +++ b/wave_train/dynamics/tdse.py @@ -19,7 +19,7 @@ class TDSE(QuantumMechanics): def __init__(self, hamilton, num_steps, step_size, sub_steps=1, solver='s2', normalize=0, - max_rank=20, repeats=20, threshold=1e-12, + max_rank=20, repeats=20, threshold=1e-12, krylov_dimension=5, save_file=None, load_file=None, compare=None): """ hamilton: instance of physical object (quantum Hamiltonian) @@ -63,6 +63,7 @@ def __init__(self, hamilton, num_steps, step_size, sub_steps=1, self.max_rank = max_rank self.repeats = repeats self.threshold = threshold + self.krylov_dimension = krylov_dimension self.save_file = save_file self.load_file = load_file self.compare = compare @@ -436,6 +437,9 @@ def update_solve(self, i): elif self.solver == 'vp': # time-dependent variational principle psi = ode.tdvp(1j*self.operator, initial_value=self.psi, step_size=self.sub_size, number_of_steps=self.sub_steps, normalize=self.normalize) + elif self.solver == 'kr': # Krylov subspace method + self.psi = ode.krylov(1j*self.operator, initial_value=self.psi, dimension=self.krylov_dimension, step_size=self.sub_size, threshold=self.threshold, + max_rank=self.max_rank, normalize=self.normalize) elif self.solver == 'qe': # Quasi-exact propagation self.psi = self.propagator @ self.psi @@ -443,7 +447,7 @@ def update_solve(self, i): sys.exit('Wrong choice of solver') # Prepare for next time step - if self.solver != 'qe': + if (self.solver != 'qe') and (self.solver != 'kr'): self.psi = psi[-1] self.psi_guess = self.psi From d06355653c7def65c8083b253824f3e9c88294a6 Mon Sep 17 00:00:00 2001 From: PGelss Date: Fri, 22 Sep 2023 11:17:51 +0200 Subject: [PATCH 19/25] included class Bath_Map_1 with test script --- test_scripts/Bath_Map_1/tdse_1.py | 97 +++++++++++++++++++ .../hamilton/{boson.py => bath_map_1.py} | 14 ++- 2 files changed, 103 insertions(+), 8 deletions(-) create mode 100644 test_scripts/Bath_Map_1/tdse_1.py rename wave_train/hamilton/{boson.py => bath_map_1.py} (86%) diff --git a/test_scripts/Bath_Map_1/tdse_1.py b/test_scripts/Bath_Map_1/tdse_1.py new file mode 100644 index 0000000..232a185 --- /dev/null +++ b/test_scripts/Bath_Map_1/tdse_1.py @@ -0,0 +1,97 @@ +from wave_train.hamilton.bath_map_1 import Bath_Map_1 +from wave_train.dynamics.tdse import TDSE +from wave_train.io.load import Load +from wave_train.io.logging import TeeLogger +from os.path import basename, splitext +import numpy as np +import matplotlib.pyplot as plt + +def exciton_tdse(batch_mode): + # Detect name of this script file (without extension) + base_name = basename(__file__) + my_file = splitext(base_name)[0] + + # logging instance: will be initialized with + # class for logging to both console and logfile + logger = None + if not batch_mode: + logger = TeeLogger(log_file=my_file + ".log") + + # Set up the excitonic Hamiltonian for a chain + hamilton = Bath_Map_1( + n_site=100, # number of sites + eta = 0.5, # coupling between TLS and first bath site + s = 1, # type of spectral density function: s<1 sub-ohmic, s=1 ohmic, s>1 super-ohmic + omega_c = 10, # cut-off frequency of spectral density function + omega_0 = 1 # eigenfrequency of the TLS + ) + + # Set up TT representation of the Hamiltonian + hamilton.get_TT( + n_basis=2, # size of electronic basis set + qtt=False # using quantized TT format + ) + + # Set up TDSE solver + dynamics = TDSE( + hamilton=hamilton, # choice of Hamiltonian, see above + num_steps=10, # number of main time steps + step_size=0.1, # size of main time steps + sub_steps=1, # number of sub steps + solver='vp', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... + normalize=0, # whether|how to normalize the solution, can be 0|2 + max_rank=5, # max rank of solution + repeats=15, # number of sweeps (implicit ODE solvers only!) + threshold=1e-8, # threshold in ALS decomposition + save_file=my_file+'.pic', # if not None, generated data will be saved to this file + load_file=None, # if not None, reference data will be loaded from this file + compare=None # How to do the comparison with reference data + ) + + + + # Set up initial state + dynamics.fundamental(list(np.eye(hamilton.n_site)[0,:])) # fundamental excitation near center of chain + # dynamics.gaussian() + # dynamics.sec_hyp(w_0=0.2) + + # Batch mode + if batch_mode: + dynamics.solve() # Solve TDSE *without* visualization + + # Interactive mode: Setup animated visualization + else: + from wave_train.graphics.factory import VisualTDSE + graphics = VisualTDSE( + dynamics=dynamics, # choice of dynamics (EoM), see above + plot_type='Populations', # select your favorite plot type + plot_expect=True, # toggle plotting of expectation values + figure_pos=(100, 50), # specifying position (x,y) of upper left of figure [in pixels] + figure_size=(1050, 450), # specifying size (w,h) of figure [in pixels] + image_file=my_file+'.png', # if not None, image (last frame) will be written to this file + movie_file=my_file+'.mp4', # if not None, animation will be written to this file + snapshots=False, # save each snapshot + frame_rate=1, # frames per second in mp4 animation file + plot_style={}, # additional plot style information + ).create() + graphics.solve() # Solve TDSE *with* visualization + + +if __name__ == '__main__': + exciton_tdse(batch_mode=False) + + + +dynamics = Load( + file_name='tdse_1', + file_type = 'pic' +) +#print(dynamics.qu_numbr) +plt.figure() +for i in range(1,100,10): + plt.plot(dynamics.qu_numbr[:,i]) +plt.show() + +#print(dynamics.qu_sig_1) + + diff --git a/wave_train/hamilton/boson.py b/wave_train/hamilton/bath_map_1.py similarity index 86% rename from wave_train/hamilton/boson.py rename to wave_train/hamilton/bath_map_1.py index 7e9b359..6776b17 100644 --- a/wave_train/hamilton/boson.py +++ b/wave_train/hamilton/bath_map_1.py @@ -3,7 +3,7 @@ from wave_train.hamilton.chain import Chain -class Boson(Chain): +class Bath_Map_1(Chain): """ Boson ... @@ -79,10 +79,8 @@ def get_2Q(self, n_dim=2): self.ground = np.zeros(n_dim) self.ground[0] = 1 self.excited = np.zeros(n_dim) - self.excited[1] = 1 - self.cat = np.zeros(n_dim) - self.cat[0] = np.sqrt(0.5) - self.cat[1] = np.sqrt(0.5) + self.excited[0] = np.sqrt(0.5) + self.excited[1] = np.sqrt(0.5) def get_SLIM(self, n_dim=2): @@ -132,10 +130,10 @@ def get_SLIM(self, n_dim=2): # define lists of ndarrays - self.S = [self.omega_0*sigmaz] + [cavity_onsite[i]*self.qu_numbr for i in range(1,self.n_site)] - self.L = [np.stack((self.eta*sigmax, self.eta*self.position), axis=-1)] + [np.stack((chain_bond[i]*self.lowering,chain_bond[i]*self.raising), axis=-1) for i in range(1,self.n_site-1)] + [None] + self.S = [-self.omega_0*sigmaz] + [cavity_onsite[i]*self.qu_numbr for i in range(1,self.n_site)] + self.L = [0.5*np.stack((self.eta*sigmax, self.eta*self.eta*sigmax), axis=-1)] + [0.5*np.stack((chain_bond[i]*self.lowering,chain_bond[i]*self.raising), axis=-1) for i in range(1,self.n_site-1)] + [None] self.I = [np.eye(3) for i in range(self.n_site)] - self.M = [None, np.stack((self.position,sigmax))] + [np.stack((self.raising,self.lowering)) for i in range(2,self.n_site)] + self.M = [None, np.stack((self.position, self.position))] + [np.stack((self.raising,self.lowering)) for i in range(2,self.n_site)] From 486845833c72255a29cd63e091fee8b56706c148 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 22 Sep 2023 16:58:40 +0200 Subject: [PATCH 20/25] Renamed Hamiltonian for exciton-phonon-coupling --- wave_train/hamilton/bath_map_1.py | 13 ++++++------- .../hamilton/{coupled.py => exc_pho_coupling.py} | 4 ++-- 2 files changed, 8 insertions(+), 9 deletions(-) rename wave_train/hamilton/{coupled.py => exc_pho_coupling.py} (97%) diff --git a/wave_train/hamilton/bath_map_1.py b/wave_train/hamilton/bath_map_1.py index 6776b17..e994908 100644 --- a/wave_train/hamilton/bath_map_1.py +++ b/wave_train/hamilton/bath_map_1.py @@ -6,10 +6,8 @@ class Bath_Map_1(Chain): """ - Boson ... - --------------------------------- - - ... + Hamiltonian for a two level system coupled to a bath + using mapping of bath degrees of freedom to a chain Parent class is chain. """ @@ -54,9 +52,10 @@ def __str__(self): omega_0 = self.omega_0 boson_str = """ ----------------------- -Hamiltonian for BOSONS ----------------------- +---------------------------------------------------- +Hamiltonian for a two level system coupled to a bath +using mapping of bath degrees of freedom to a chain +---------------------------------------------------- Coupling between TLS and first bath site (eta) = {} Type of spectral density function (s) = {} diff --git a/wave_train/hamilton/coupled.py b/wave_train/hamilton/exc_pho_coupling.py similarity index 97% rename from wave_train/hamilton/coupled.py rename to wave_train/hamilton/exc_pho_coupling.py index f64db79..be4f31c 100644 --- a/wave_train/hamilton/coupled.py +++ b/wave_train/hamilton/exc_pho_coupling.py @@ -8,7 +8,7 @@ class Coupled(Chain): """ - Coupled electron-phonon dynamics + exciton-phonon coupling dynamics -------------------------------- for a chain of exciton sites, connected by harmonic oscillators @@ -169,7 +169,7 @@ def __str__(self): coupled_str = """ ----------------------------- -Coupling EXCITONS and PHONONS +EXCITONS - PHONONS - Coupling ----------------------------- Linear tuning constant (chi), localized = {} From ab69151942289e16ea45dea1c2f0c446a37c0a54 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 22 Sep 2023 17:49:36 +0200 Subject: [PATCH 21/25] Adapted test scripts for exciton-phonon-coupling to latest code changes --- test_scripts/Bath_Map_1/tdse_1.py | 25 ++++++++----------- .../{Coupled => Exc_Pho_Coupling}/qcmd_1.py | 4 +-- .../{Coupled => Exc_Pho_Coupling}/qcmd_2.py | 4 +-- .../{Coupled => Exc_Pho_Coupling}/qcmd_3.py | 4 +-- .../{Coupled => Exc_Pho_Coupling}/qcmd_4.py | 4 +-- .../{Coupled => Exc_Pho_Coupling}/tdse_1.py | 2 +- .../tdse_qe_1.py | 6 ++--- .../tdse_qe_2.py | 6 ++--- .../tdse_qe_3.py | 4 +-- .../{Coupled => Exc_Pho_Coupling}/tise_1.py | 4 +-- wave_train/hamilton/bath_map_1.py | 3 ++- wave_train/hamilton/exc_pho_coupling.py | 12 ++++----- 12 files changed, 38 insertions(+), 40 deletions(-) rename test_scripts/{Coupled => Exc_Pho_Coupling}/qcmd_1.py (95%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/qcmd_2.py (95%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/qcmd_3.py (95%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/qcmd_4.py (95%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/tdse_1.py (96%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/tdse_qe_1.py (93%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/tdse_qe_2.py (94%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/tdse_qe_3.py (95%) rename test_scripts/{Coupled => Exc_Pho_Coupling}/tise_1.py (95%) diff --git a/test_scripts/Bath_Map_1/tdse_1.py b/test_scripts/Bath_Map_1/tdse_1.py index 232a185..f1d4786 100644 --- a/test_scripts/Bath_Map_1/tdse_1.py +++ b/test_scripts/Bath_Map_1/tdse_1.py @@ -6,7 +6,7 @@ import numpy as np import matplotlib.pyplot as plt -def exciton_tdse(batch_mode): +def bath_map_1_tdse(batch_mode): # Detect name of this script file (without extension) base_name = basename(__file__) my_file = splitext(base_name)[0] @@ -19,11 +19,11 @@ def exciton_tdse(batch_mode): # Set up the excitonic Hamiltonian for a chain hamilton = Bath_Map_1( - n_site=100, # number of sites + n_site=50, # number of sites eta = 0.5, # coupling between TLS and first bath site s = 1, # type of spectral density function: s<1 sub-ohmic, s=1 ohmic, s>1 super-ohmic omega_c = 10, # cut-off frequency of spectral density function - omega_0 = 1 # eigenfrequency of the TLS + omega_0 = 1 # eigenfrequency of the TLS ) # Set up TT representation of the Hamiltonian @@ -36,19 +36,17 @@ def exciton_tdse(batch_mode): dynamics = TDSE( hamilton=hamilton, # choice of Hamiltonian, see above num_steps=10, # number of main time steps - step_size=0.1, # size of main time steps + step_size=0.1, # size of main time steps sub_steps=1, # number of sub steps solver='vp', # can be 'se' (symmetrized Euler) or 'sm' (Strang-Marchuk splitting) or ... normalize=0, # whether|how to normalize the solution, can be 0|2 max_rank=5, # max rank of solution repeats=15, # number of sweeps (implicit ODE solvers only!) - threshold=1e-8, # threshold in ALS decomposition + threshold=1e-8, # threshold in ALS decomposition save_file=my_file+'.pic', # if not None, generated data will be saved to this file load_file=None, # if not None, reference data will be loaded from this file compare=None # How to do the comparison with reference data ) - - # Set up initial state dynamics.fundamental(list(np.eye(hamilton.n_site)[0,:])) # fundamental excitation near center of chain @@ -70,7 +68,7 @@ def exciton_tdse(batch_mode): figure_size=(1050, 450), # specifying size (w,h) of figure [in pixels] image_file=my_file+'.png', # if not None, image (last frame) will be written to this file movie_file=my_file+'.mp4', # if not None, animation will be written to this file - snapshots=False, # save each snapshot + snapshots=False, # save each snapshot frame_rate=1, # frames per second in mp4 animation file plot_style={}, # additional plot style information ).create() @@ -78,10 +76,9 @@ def exciton_tdse(batch_mode): if __name__ == '__main__': - exciton_tdse(batch_mode=False) - - - + bath_map_1_tdse(batch_mode=False) + + dynamics = Load( file_name='tdse_1', file_type = 'pic' @@ -93,5 +90,5 @@ def exciton_tdse(batch_mode): plt.show() #print(dynamics.qu_sig_1) - - + + diff --git a/test_scripts/Coupled/qcmd_1.py b/test_scripts/Exc_Pho_Coupling/qcmd_1.py similarity index 95% rename from test_scripts/Coupled/qcmd_1.py rename to test_scripts/Exc_Pho_Coupling/qcmd_1.py index 773111e..dc024f6 100644 --- a/test_scripts/Coupled/qcmd_1.py +++ b/test_scripts/Exc_Pho_Coupling/qcmd_1.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.qcmd import QCMD from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_qcmd(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=15, # number of sites periodic=False, # periodic boundary conditions homogen=True, # homogeneous chain/ring diff --git a/test_scripts/Coupled/qcmd_2.py b/test_scripts/Exc_Pho_Coupling/qcmd_2.py similarity index 95% rename from test_scripts/Coupled/qcmd_2.py rename to test_scripts/Exc_Pho_Coupling/qcmd_2.py index b998ac9..54a6026 100644 --- a/test_scripts/Coupled/qcmd_2.py +++ b/test_scripts/Exc_Pho_Coupling/qcmd_2.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.qcmd import QCMD from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -17,7 +17,7 @@ def coupled_qcmd(batch_mode): # Define properties of chain/ring system # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=15, # number of sites periodic=False, # periodic boundary conditions homogen=True, # Homogeneous chain/ring diff --git a/test_scripts/Coupled/qcmd_3.py b/test_scripts/Exc_Pho_Coupling/qcmd_3.py similarity index 95% rename from test_scripts/Coupled/qcmd_3.py rename to test_scripts/Exc_Pho_Coupling/qcmd_3.py index 5dbeae1..6f7272c 100644 --- a/test_scripts/Coupled/qcmd_3.py +++ b/test_scripts/Exc_Pho_Coupling/qcmd_3.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.qcmd import QCMD from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_qcmd(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=41, # number of sites periodic=True, # periodic boundary conditions homogen=True, # Homogeneous chain/ring diff --git a/test_scripts/Coupled/qcmd_4.py b/test_scripts/Exc_Pho_Coupling/qcmd_4.py similarity index 95% rename from test_scripts/Coupled/qcmd_4.py rename to test_scripts/Exc_Pho_Coupling/qcmd_4.py index b1bfd16..add6bda 100644 --- a/test_scripts/Coupled/qcmd_4.py +++ b/test_scripts/Exc_Pho_Coupling/qcmd_4.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.qcmd import QCMD from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_qcmd(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=41, # number of sites periodic=False, # periodic boundary conditions homogen=True, # Homogeneous chain/ring diff --git a/test_scripts/Coupled/tdse_1.py b/test_scripts/Exc_Pho_Coupling/tdse_1.py similarity index 96% rename from test_scripts/Coupled/tdse_1.py rename to test_scripts/Exc_Pho_Coupling/tdse_1.py index 9a43050..cb6fab0 100644 --- a/test_scripts/Coupled/tdse_1.py +++ b/test_scripts/Exc_Pho_Coupling/tdse_1.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.tdse import TDSE from wave_train.io.logging import TeeLogger from os.path import basename, splitext diff --git a/test_scripts/Coupled/tdse_qe_1.py b/test_scripts/Exc_Pho_Coupling/tdse_qe_1.py similarity index 93% rename from test_scripts/Coupled/tdse_qe_1.py rename to test_scripts/Exc_Pho_Coupling/tdse_qe_1.py index 88685d2..56a8d95 100644 --- a/test_scripts/Coupled/tdse_qe_1.py +++ b/test_scripts/Exc_Pho_Coupling/tdse_qe_1.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.tdse import TDSE from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_qe_1(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=3, # number of sites periodic=False, # periodic boundary conditions homogen=True, # homogeneous chain/ring @@ -39,7 +39,7 @@ def coupled_qe_1(batch_mode): # Set up TDSE solver dynamics = TDSE( - hamilton=hamilton, # choice of Hamiltonian, see above + hamilton=hamilton, # choice of Hamiltonian, see above num_steps=50, # number of main time steps step_size=20, # size of main time steps sub_steps=None, # dummy diff --git a/test_scripts/Coupled/tdse_qe_2.py b/test_scripts/Exc_Pho_Coupling/tdse_qe_2.py similarity index 94% rename from test_scripts/Coupled/tdse_qe_2.py rename to test_scripts/Exc_Pho_Coupling/tdse_qe_2.py index 32ce4bb..157375f 100644 --- a/test_scripts/Coupled/tdse_qe_2.py +++ b/test_scripts/Exc_Pho_Coupling/tdse_qe_2.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.tdse import TDSE from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_qe_2(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=3, # number of sites periodic=False, # periodic boundary conditions homogen=True, # homogeneous chain/ring @@ -39,7 +39,7 @@ def coupled_qe_2(batch_mode): # Set up TDSE solver dynamics = TDSE( - hamilton=hamilton, # choice of Hamiltonian, see above + hamilton=hamilton, # choice of Hamiltonian, see above num_steps=50, # number of main time steps step_size=20, # size of main time steps sub_steps=20, # number of sub steps diff --git a/test_scripts/Coupled/tdse_qe_3.py b/test_scripts/Exc_Pho_Coupling/tdse_qe_3.py similarity index 95% rename from test_scripts/Coupled/tdse_qe_3.py rename to test_scripts/Exc_Pho_Coupling/tdse_qe_3.py index 1412ee3..08d6102 100644 --- a/test_scripts/Coupled/tdse_qe_3.py +++ b/test_scripts/Exc_Pho_Coupling/tdse_qe_3.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.qcmd import QCMD from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_qe_3(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=3, # number of sites periodic=False, # periodic boundary conditions homogen=True, # homogeneous chain/ring diff --git a/test_scripts/Coupled/tise_1.py b/test_scripts/Exc_Pho_Coupling/tise_1.py similarity index 95% rename from test_scripts/Coupled/tise_1.py rename to test_scripts/Exc_Pho_Coupling/tise_1.py index b1b22ba..3160fd3 100644 --- a/test_scripts/Coupled/tise_1.py +++ b/test_scripts/Exc_Pho_Coupling/tise_1.py @@ -1,4 +1,4 @@ -from wave_train.hamilton.coupled import Coupled +from wave_train.hamilton.exc_pho_coupling import Exc_Pho_Coupling from wave_train.dynamics.tise import TISE from wave_train.io.logging import TeeLogger from os.path import basename, splitext @@ -15,7 +15,7 @@ def coupled_tise(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled excitonic-phononic Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=5, # number of sites periodic=True, # periodic boundary conditions homogen=True, # homogeneous chain/ring diff --git a/wave_train/hamilton/bath_map_1.py b/wave_train/hamilton/bath_map_1.py index e994908..fede487 100644 --- a/wave_train/hamilton/bath_map_1.py +++ b/wave_train/hamilton/bath_map_1.py @@ -74,7 +74,8 @@ def get_2Q(self, n_dim=2): # Call corresponding method from superclass Chain super().get_2Q(n_dim) - # Set up ground and first excited state + # Set up ground and excited state + # Assuming that the two-level-system is in a Schroedinger CAT state self.ground = np.zeros(n_dim) self.ground[0] = 1 self.excited = np.zeros(n_dim) diff --git a/wave_train/hamilton/exc_pho_coupling.py b/wave_train/hamilton/exc_pho_coupling.py index be4f31c..c1244ad 100644 --- a/wave_train/hamilton/exc_pho_coupling.py +++ b/wave_train/hamilton/exc_pho_coupling.py @@ -6,10 +6,10 @@ from wave_train.hamilton.exciton import Exciton -class Coupled(Chain): +class Exc_Pho_Coupling(Chain): """ - exciton-phonon coupling dynamics - -------------------------------- + exciton-phonon coupling + ----------------------- for a chain of exciton sites, connected by harmonic oscillators optionally with periodic boundaries @@ -168,9 +168,9 @@ def __str__(self): tau = self.tau coupled_str = """ ------------------------------ -EXCITONS - PHONONS - Coupling ------------------------------ +--------------------------- +EXCITON - PHONON - Coupling +--------------------------- Linear tuning constant (chi), localized = {} Linear tuning constant (rho), non-symmetric = {} From 63848a9d04e6c31183bbf8ed25727fa92439bf6b Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Fri, 22 Sep 2023 18:45:21 +0200 Subject: [PATCH 22/25] Encoding the Krylov dimension in the name of the integrator --- wave_train/dynamics/tdse.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/wave_train/dynamics/tdse.py b/wave_train/dynamics/tdse.py index 26234c9..8752cda 100644 --- a/wave_train/dynamics/tdse.py +++ b/wave_train/dynamics/tdse.py @@ -19,7 +19,7 @@ class TDSE(QuantumMechanics): def __init__(self, hamilton, num_steps, step_size, sub_steps=1, solver='s2', normalize=0, - max_rank=20, repeats=20, threshold=1e-12, krylov_dimension=5, + max_rank=20, repeats=20, threshold=1e-12, save_file=None, load_file=None, compare=None): """ hamilton: instance of physical object (quantum Hamiltonian) @@ -63,7 +63,6 @@ def __init__(self, hamilton, num_steps, step_size, sub_steps=1, self.max_rank = max_rank self.repeats = repeats self.threshold = threshold - self.krylov_dimension = krylov_dimension self.save_file = save_file self.load_file = load_file self.compare = compare @@ -437,8 +436,8 @@ def update_solve(self, i): elif self.solver == 'vp': # time-dependent variational principle psi = ode.tdvp(1j*self.operator, initial_value=self.psi, step_size=self.sub_size, number_of_steps=self.sub_steps, normalize=self.normalize) - elif self.solver == 'kr': # Krylov subspace method - self.psi = ode.krylov(1j*self.operator, initial_value=self.psi, dimension=self.krylov_dimension, step_size=self.sub_size, threshold=self.threshold, + elif self.solver in ['k2', 'k4', 'k6', 'k8']: # Krylov subspace method + self.psi = ode.krylov(1j*self.operator, initial_value=self.psi, dimension=int(self.solver[1]), step_size=self.sub_size, threshold=self.threshold, max_rank=self.max_rank, normalize=self.normalize) elif self.solver == 'qe': # Quasi-exact propagation self.psi = self.propagator @ self.psi @@ -447,7 +446,7 @@ def update_solve(self, i): sys.exit('Wrong choice of solver') # Prepare for next time step - if (self.solver != 'qe') and (self.solver != 'kr'): + if self.solver not in ['qe', 'k2', 'k4', 'k6', 'k8']: self.psi = psi[-1] self.psi_guess = self.psi From c85ccd1cae9e1eceafb3ce0b3bf2cb635b455c7b Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Mon, 25 Sep 2023 15:49:29 +0200 Subject: [PATCH 23/25] Added loop over sub-steps for Krylov propagations --- wave_train/dynamics/tdse.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/wave_train/dynamics/tdse.py b/wave_train/dynamics/tdse.py index 8752cda..536534d 100644 --- a/wave_train/dynamics/tdse.py +++ b/wave_train/dynamics/tdse.py @@ -437,8 +437,9 @@ def update_solve(self, i): psi = ode.tdvp(1j*self.operator, initial_value=self.psi, step_size=self.sub_size, number_of_steps=self.sub_steps, normalize=self.normalize) elif self.solver in ['k2', 'k4', 'k6', 'k8']: # Krylov subspace method - self.psi = ode.krylov(1j*self.operator, initial_value=self.psi, dimension=int(self.solver[1]), step_size=self.sub_size, threshold=self.threshold, - max_rank=self.max_rank, normalize=self.normalize) + for step in range(self.sub_steps): + self.psi = ode.krylov(1j*self.operator, initial_value=self.psi, dimension=int(self.solver[1]), step_size=self.sub_size, threshold=self.threshold, + max_rank=self.max_rank, normalize=self.normalize) elif self.solver == 'qe': # Quasi-exact propagation self.psi = self.propagator @ self.psi From 15fc3c9b008575655d136735be16076d203f0690 Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Tue, 26 Sep 2023 15:23:38 +0200 Subject: [PATCH 24/25] After the merging, this is intended to become version 1.0.11 --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 4f57566..feb939d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = wave_train -version = 1.0.10 +version = 1.0.11 author = Jerome Riedel, Patrick Gelß, Burkhard Schmidt author_email = burkhard.schmidt@fu-berlin.de description = Numerical quantum mechanics of chain-like systems based on tensor trains From 5332665d4978b3182a324ff2ab19d0c14628fbfc Mon Sep 17 00:00:00 2001 From: Burkhard Schmidt Date: Wed, 27 Sep 2023 16:02:11 +0200 Subject: [PATCH 25/25] Minor fix: hamilton = Exc_Pho_Coupling --- test_scripts/Exc_Pho_Coupling/tdse_1.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_scripts/Exc_Pho_Coupling/tdse_1.py b/test_scripts/Exc_Pho_Coupling/tdse_1.py index cb6fab0..b8464d6 100644 --- a/test_scripts/Exc_Pho_Coupling/tdse_1.py +++ b/test_scripts/Exc_Pho_Coupling/tdse_1.py @@ -15,7 +15,7 @@ def coupled_tdse(batch_mode): logger = TeeLogger(log_file=my_file + ".log") # Set up the coupled exciton-phonon Hamiltonian for a chain - hamilton = Coupled( + hamilton = Exc_Pho_Coupling( n_site=15, # number of sites periodic=False, # periodic boundary conditions homogen=True, # homogeneous chain/ring