tmancal74 / quantarhei

Open Quantum System Theory for Molecular Systems

QUANTArhei: Open Quantum System Theory for Molecular Systems

Quantarhei is a Molecular Open Quantum Systems Simulator written predominantly in Python. Its name is derived from the famous aphorism "Panta rhei" of the Greek philosopher Heraclitus of Ephesus. "Panta rhei" means "Everything flows" or "Everything is in flux" which is quite fitting when you change Panta into Quanta.

In "Quantarhei" the last four letters ("rhei") should be written in Greek, i.e. (using LateX convention) "\rho \epsilon \iota".

---

Quantarhei is in flux, but it already provides helper classes to define molecules, their aggregates and their interaction with external environment. It can calculate absorption spectra of individual molecules and their aggregates and excitation energy transfer dynamics using various types of Redfield and Foerster theories.

Quantarhei provides Python code (optimized with Numpy) for all its implemented methods and theories, and allows extensions and replacements of the reference Python code with optimised routines written in C, Fortran or other lower level languages.

In the first development stage, we concentrate on bringing to you tools to quickly build essential components of a quantum mechanical simulation, such as Hamiltonian, relaxation tensors, various initial conditions for density matrix etc.

Quantarhei is at its experimental stage. Current version is

Quantarhei is available in source form on GitHub and from PyPI for installation with the pip command.

Acknowledgements

The work on Quantarhei is supported by

Neuron Fund for Support of Science

through the Impuls grant in physics 2014 (2015-2017)

and

Czech Science Foundation (GACR)

through grants: 14-25752S (2014-2016), 17-22160S (2017-2019) and 18-18022S (2018-2020)

New in 0.0.66

for users:

• Quantarhei starts enforcing usage of 'energy_units' context with some function. When energy values are set, it will be compulsary to use the 'energy_units' context manager. The introduction of this feature may take several Quantarhei versions.
• bug fixes

New in 0.0.65

for users:

• bug fixes

New in 0.0.64

for users:

• use cases introduced for input files
• bug fixes

New in 0.0.63

for users:

• bug fixes

New in 0.0.62

for users:

• _input_file_ variable defined in scripts run through qrhei run command
• assert_version() function defined
• bug fixes

New in 0.0.61

for users:

• bug fixes

New in 0.0.60

for users:

• bug fixes

New in 0.0.59

for users:

• bug fixes

New in 0.0.58

for users:

• Many bug fixes

New in 0.0.57

for users:

• Bug fixes

New in 0.0.56

for users:

• Bug fixes

New in 0.0.55

for users:

• A bug in 2D lineshapes of ESA calculated by MockTwoDSpectrumCalculator fixed
• Other bug fixes

New in 0.0.54

for users:

• TwoDSpectra can be added if they are of the same type
• Bug fixes

New in 0.0.53

for users:

• Improved handling of rotating wave approximation (RWA) for density matrix and the state vector
• Bug fixes

New in 0.0.52

for users:

• Improved control over parallelization from command line by qrhei options
• Parallelization over multiple nods
• Bug fixes

New in 0.0.51

for users:

• Bug fixes

New in 0.0.50

for users:

• Some improvement of the qrhei driver
• Improved logging capabilities; standard print function can be replaced printlog function
• Bug fixes

New in 0.0.49

for users:

• Runnable .yaml configuration files
• Better diagnostics of errors occuring while running a script by qrhei driver
• New qrhei subcommand file, which shows information about files saved by Quantarhei
• Bug fixes

for developers:

• Started work on a script compiler qtask

New in 0.0.48

for users:

• Bug fixes

New in 0.0.47

for users:

• Bug fixes
• 2D spectra can be shifted in its axis by less then the frequency step (interpolated shift)
• New helper class Input to simplify configuration of user scripts by ".yaml" or ".json" files

New in 0.0.46

for users:

• A bug introduced in 0.0.45 in 2D spectrum calculations now fixed
• Most of the classes can now be copied (.copy() for a shallow copy, .deepcopy() for a deep, recursive copy)
• Improvements of 2D calculations with dressed states
• Vibronic example of 2D spectrum calculation extended
• Class migration: TwoDSpectrum -> TwoDResponse; TwoDSpectrumContainer -> TwoDResponseContainer
• New classes TwoDSpectrum and TwoDSpectrumContainer are simpler and contain only one type of spectra
• Class migration: MockTwoDSpectrum -> Mocj
• Definitions of positive and negative frequencies in omega2-frequency maps changed to correspond to literature

For developers:

• many constants describing non-linear response, 2D spectra and data are defined on the highest quantarhei import level
• All Saveable objects now have a convenience methods copy() for shallow copying, deepcopy() for deep copying and scopy() for a deep copy by saving to and loading the object from a temporary file

New in 0.0.45

For users:

• Improved ability to calculate and analyze 2D frequency maps
• get_Fourier_transform method of DFunction accepts a windowing function, which works the same as the one of TwoDSpectrumContainer
• Calculation of effective lineshape pump-probe spectra
• Saving of TwoDSpectrum data into text files (.dat, .txt), numpy formats (.npy, .npz) and Matlab files (.mat) is enabled
• Problem which caused Redfield and Lindblad operators not to transform to correct basis when represented by operators was fixed
• Problem which caused Liouville pathways to be calculated with site basis evolution superoperator was fixed
• operator_factory class of Harmonic oscillator now fixed to return correct shift operator for complex inputs
• operator_factory is now available from quantarhei.models.HarmonicOscillator package
• Some small bug fixes

For developers:

• Failing documentation compilation fixed and adjusted to new version of matplotlib

New in 0.0.44

For users:

• Basic implementation of HEOM
• Some bug fixes

New in 0.0.43

For users:

• PureDephasing super-operator to allow additional pure dephasing for realistic lineshapes in effective lineshape description of time-resolved experiments
• Empty relaxation superopetator (as an empty Lindblad form) introduced (as a temporary fix to allow pure dephasing dynamics only)
• Consistent calculation of pure dephasing of non-optical coherence elements of the density matrix from effective lineshape theory (including electronic only dephasing in vibrational-electronic systems)
• Some bug fixes

New in 0.0.42

For users:

• Improved effective lineshapes for 2D spectrum calculations
• Calculation of absorption spectrum using first order Liouville pathways
• Some bug fixes including an frequency factor in absorption spectrum

New in 0.0.41

For users:

• Some bug fixes
• Better Louville pathway manipulation features

New in 0.0.40

For users:

• Some bug fixes
• Minor new features

New in 0.0.39

For users:

• Some bug fixes

New in 0.0.38

For users:

• Some bug fixes

New in 0.0.37

For users:

• Some bug fixes

For developers

• Some unused files removed
• More precise dependencies on other packages specified in setup

New in 0.0.36

For users:

• Quantarhei now available also as a conda package
• Recommended installation procedure documented
• TwoDSpectrum class revised - new method names, better storage model (keeps track of rephasing and non-rephasing part, groups of pathways associated with different processes when required, stores different pathways separately when required)
• Improved TwoDSpectrumContainer (can hold a group of spectra identified by an arbitrary ValueAxis (most notably TimeAxis and FrequencyAxis), integer index or list of strings). Copies the new storage improvement on TwoDSpectrum.
• labsetup class changed to LabSetup and extended by information about pulse profiles and spectra. labsetup is left as deprecated for compatibility
• Fourier transform of 2D spectra in t2, via TwoDSpectrumContainer; also enables FFT with window function
• Functions of ValueAxis introduced in a special module; Tukey window function for FFT in waiting time is one of them
• SuperOperator is BasisManaged; basis management is solved for both time-dependent and time-independent super operators
• RelaxationTensor now inherits from SuperOperator and it is BasisManaged through that inheritance
• EvolutionSuperOperator tested, documented and it is BasisManaged
• EvolutionSuperOperator’s method apply() can be applied with time argument which is of type TimeAxis type, float or array of floats; returns DensityMatrix or DensityMatrixEvolution
• Quantarhei driver qrhei changes format: use ‘qrhei run scriptname’ to run scripts and consult the -h option of ‘qrhei run’; parallel runs untested in this version
• Documentation contains a description of the concept of “user”, “advanced”, and “expert” levels of classes in Quantarhei.
• List of classes completely covered by documentation and doctests included in on-line documentation
• Classes Mode, SubMode, Molecule, TwoDSpectrumContainer completely documented
• Documentation enhanced
• Countless small improvements and bug fixes

For developers:

• Code of conduct file now in the root directory of the package
• Absorption spectroscopy related classes now organized in one file per class fashion so that automatic documentation is easier to read
• New subpackage quantarhei.testing united all custom functions that support testing. It includes feature.py module previously found in quantarhei.dev subpacked (now removed) and a behave.py module which supports tests with behave package
• Behave package is now used for some tests (in particular for tests of the “qrhei” driver). Future acceptance tests should preferentially be written with this package
• New helper script “ghenerate” autogenerates Python step files for tests with ‘behave’ package from the Gherkin feature files

New in 0.0.35

For users:

• Method get_DensityMatrix() of the Aggregate class improved. It accepts some new options which makes specification of desired density matrix more flexible
• Experimental implementation of circular and linear dichroisms and fluorescence spectra
• Documentation is now available on readthedocs.org. A badge which informations about the status of automatic documentation builds was added to README
• Many small improvements and bug fixes

For developers:

• The code is now hosted on travis-ci.com and the builds are tested after every commit. Corresponding badge has been added to README
• The code is now hosted on codecov.com and its coverage by tests is measured. Corresponding badge showing the coverage has beed added to README

New in 0.0.34

For users

• Some issues with addition of bath correlation functions was fixed
• First entry in a database of literature bath correlation functions was created: the vibrational part of the FMO spectral density from Wendling et al., (2004)
• Aggregate can return a matrix of Franck-Condon factors (get_FC_factor_matrix())
• Aggregate can transform excited state site-basis shifted vibrational representation of an arbitrary operator to the unshifted (ground state) one (transform_2_unshifted(A, inverse=True/False) )
• Several new tested examples
• RelaxationTensors (Redfield, Foerster, Lindblad, etc.) can now be multiplied by a constant or added (addition only if they are in tensor, i. e. not in operator, form)
• Tested examples can be fetched into IPython notebook or Python/IPython console by %example magic command or fetch_example function from quantarhei.wizard.magic module
• Small improvements and bug fixes

New in 0.0.33

For users:

• Evolution superoperators for relaxation tensors with constant coefficients (EvolutionSuperOperator class)
• Liouville pathway analysis including relaxation pathways (in Aggregate class)
• Small improvements and bug fixes

For developers:

• Aggregate class is broken into smaller pieces which snowball the functionality. Basic class is AggregateBase; new functions of this powerful class are defined in separate child classes. Aggregate class inherits from the whole chain of classes
• quantarhei.REAL and quantarhei.COMPLEX types should be now used for numpy arrays throughout the package. These types can be controlled and with it the used numerical precision and memory needs

New in 0.0.32

For users:

• Electronic Lindblad form for vibronic Frenkel exciton model
• Propagation with relaxation tensor (in particular Redfield and Time-dependent Redfield) in operator representation (where applicable it is much faster than with the tensorial representation)
• Redfield tensor and Time-dependent Redfield tensor can be calculated for a model with arbitrary number of vibrational states
• Aggregate can vibrationally trace arbitrary operator defined on its Hilbert space
• Small improvements and bug fixes

New in version 0.0.31

For users:

• Arbitrary time independent Lindblad form
• quantarhei.wizard module which contains IPython magic commands and some helpful Python console commands
• Simulation templates which can be fetched into IPython notebooks or console by %template magic command (IPython) or fetch_template (console and IPython)
• Part of the test suit available for installed Quantarhei package
• Some small improvements and bug fixes

For developers:

• Makefile is back in the package root directory
• examples directory depleted in favor of quantarhei/wizard/examples directory
• New tests under quantarhei/tests directory (mostly unit tests which contain plots)
• pytest required to run newtests with matplotlib plots

quantarhei's Issues

Reimplement absorption spectrum using the DFunction and its Fourier transform capabilities.

AbsSpect class has to be more based on DFunction and its Fourier transform capabilities. It also should use the symmetry inbuilt into the definition of the upper-half type of the TimeAxis.

Because of a fix in Correlation function, absorption tests now fail

Test data of absorption has to be recreated. They contained a wrong sign of the imaginary part of the correlation function.

Aggregate objects have to be able to create a default SystemBathInteraction object from the available molecules

If implemented, this enhancement with greatly simplify the definition of a problem where uncorrelated energy gap fluctuations are assumed.

Units management in FrequencyAxis is not complete

Some components of the FrequencyAxis objects do not react to the management of units. This causes confusion when calculating FT and similar stuff.

Compare Redfield population rates with Redfield tensor calculations (for FG and RC projects)

Redfield rates are tested against analytical results. Redfield tensor should be tested against the rates.

Fix absorption spectrum calculation for Molecules and Aggregates (for FG project)

Absorption spectrum for an aggregate is calculated in a wrong way. There is some shift in the indices of states which leads to a wrong assignment of dipole moments to excitonic states

Enable trimming of 2D spectra to save memory and diskspace

Calculated 2D spectra are too large for storage and they contain mostly many zeros on the edges. Provide function trim_to(window) for the TwoDSpectrum object. Something similar should be possible also for the TwoDSpectrumContainer

PopulationPropagator and ReducedDensityMatrixPropagator do not give consistent time evolutions

Implement Fluorescence and Photon Counting (for RC project)

Existing project needs only relaxation rates for calculation of time-dependent fluorescence. It is time to implement it.

Fix plotting of DFunction so that it can be used to produce movies (see example for pump-probe)

Unify the plotting options for DFunction (AbsSpectrum, PumpProbeSpectrum etc.) and TwoDSpectrum. Make shore that one can always make make movies like for 2D spectra.

Implement correlation function and spectra density definitions.

• Correlation function is not aware of its dimension of energy^2
• SpectralDensity class is not synchronised with CorrelationFunction class (as it should according to wiki)
• Related: DFunction interpolation does not work for complex functions! Only real part is interpolated.

Core.frequency for odd numer

if (self.length % 2) != 0:
raise Exception("Cannot create upper-half TimeAxis"
+ " from an odd number of points")
I feel uncomfortable
Try to compare with

https://github.com/whaley-group-berkeley/QSpectra/blob/master/qspectra/simulate/utils.py
Also for symmetry
FIXME: No choice of symmetry provided !!!!
for k in range(0, t.length-1):
yy[w.length-k-1] = numpy.conj(y[k+1])
While in Qspectra just add zeros

Retest the AbsSpect class

Test the results visually and create a new test suite.

Cyclic references of the objects to themselves have to be handled

Currently, cyclic references throw an Exception. They should rather be correctly saved and restored upon loading

Move fitting into DFunction

Fitting of absorption spectra should be provided by DFunction class and AbsSpectrum class

Put AbsSpect under units management

AbsSpect class has several methods which need to be put under units management.

• calculate(rwa) needs to understand rwa in current units
• plotting needs to react to frequency units

Setting of relaxation in vibrational modes

When a mode is defined for a molecule, its interaction with the environment has to be specified explicitly for each electronic state of the molecule. This is too complicated. Most often the interaction will the the same for all states. An update should complete the setting automatically for all electronic states.

In some tests the units management works only partially.

It seems that in the acceptance tests of creation of Redfield tensor (run by Paver) when printing an energy managed object, such as Hamiltonian, the values are not properly converted to current units. Units conversion at definition works.

Implement basis control for relaxation tensors

This is tricky, because relaxation tensors are often created in a very specific basis. Straightforward implementation should be easy, but there must be some special usage of the basis control!

Enforcement of units and basis context usage

It is desirable not to allow functions that set physical values to be used outside the units context (energy_units etc.). These functions should raise an Exception if used without the context. Similarly, some functions should not be used inside the eigenbasis_of context. Again, an Exception should be raised if this context is applied on the function.

Configuration file qrconf.py is not respected in Jupyter notebooks

It seems that in Jupyter notebook, the configuration file is aways rewritten by default, and new settings do not affect the execution.

Convert AbsSpectrum.save to use hdf5 format

AbsSpectrumContainer already used hdf5; use it as default

Implement measure_reorganization_energy method in SpectralDensity and CorrelationFunction classes

measure_reorganization_energy() method is important for tracking the consistency of the two classes, especially if some of the quantities are calculated by numerical transforms.

Make CorrelationFunction.get_SpectralDensity() return a numerically calculated spectral density

a Numerical rather than a "per definition" version of the spectral density should be returned by this method. The "per definition" version is obtained by constructing SpectralDensity object directly and sometimes it can be obtained analytically. Comparing the numerical and analytical form is a good check.

Implement secularization through inheritance from Secular class

Issues #75 and #65 can be fixed by implementing a class (Secular) which provides basic framework for secularization of an arbitrary SuperOperator included Relaxation tensors. In case of secularization, the SuperOperator will keep a matrix of relaxation rates and matrix of dephasing rates (in case of time-dependent operators, it will keep time evolution of the rates) and classes which handle propagation will use this information to efficiently propagate. Also EvolutionSuperOperator can be secularized or calculated in secular approximation

Affected classes:

SuperOperator
EvolutionSuperOperator
RelaxationTensor (and its derivatives)

Secular approximation not available for relaxation tensor in operator form

Operator form of relaxation operators is now default, but secularization fails in this case. Secular version of relaxation operator should be implemented also for operator form.

Reaching correct correlation function for a given state from SystemBathInteraction object does not work

The get_coft(n,m) method of the SBI returns wrong correlation function. Overall, there has to be made some simple logic of getting correlation functions from CorrelationFunctionMatrix and SBI, and it has to be described somewhere!

Tutorials

I wish tutorial to cover textbook
Molecular Excitation Dynamics
and Relaxation

This will be so helpful for beginners

Time dependent versions of Redfield, Foerster and their combination.

Complete the control of relaxation theories from Aggregate object.

Some correlation function and spectral density types do not go through automatic units conversion of their parameters

CorrelationFunction and SpectralDensity objects need to be rechecked for automated conversion of the parameters of energy dimension. So types of predefined functions do not do this conversion properly.

Rotating Wave Approximation in propagators has to be used transparently

Current implementation of the propagators does not give use any control over the use of RWA. We need to implement switches to allow the user to prevent RWA, or enable output of the propagated density matrix in RWA. The default behavior should be that RWA is used if Hamiltonian contains information about the right RWA frequencies, but the calculated density matrix should be outputted including the fast frequencies.

Create acceptance tests for absorption spectrum calculations

There are no acceptance tests for AbsSpect class. Since now the code is tested manually and seems to work, it would be the right time to "freeze" some of its results into a file against which one could compare in future (when a more general method of calculation of the spectrum is implemented).

OpenSystem does not know how to create Kubo-Tanimura hierarchy

Kubo-Tanimura hierarchy and the correspoding propagator can be created from the configured SBI and Hamiltonian. OpenSystem object has these attributes and it therefore deserves an implementation of the methods, such as get_KTHierarchy and get_KTHierarchyPropagator to enable an easy application of the HEOM theory on the Molecule and Aggregate objects.

Basis control on Relaxation tensor requires some more elaboration. Currently, there are too many basis changes enforced!

Because Redfield tensor HAS TO be calculated in the excitonic basis, and the Hamiltonian HAS TO be submitted in the local basis, we end up with unnecessary transformations of the whole tensor right after its creation. This has to be circumvented! Suggestion is to implement basis management context managers which KNOW about excitonic and diabatic basis and avoid unnecessary transformation enforce by the fact that current basis context managers blindly transform when basis is changed.

Evolution superoperator is not basis aware

Evolution superoperator calculation does not react to the eigenbasis_of context. It should be make aware of it, and it should be used in diagonal basis when applied in calculation of Liouville pathways.

Implement a global set_current_units(dict) method of the Manager class

Implement

Manager().set_current_units(dict)

method, mentioned in the wiki Management of Units section

Set "as_operators" default for Relaxation tensors to True

If this default setting is changed for Redfield tensor (and its derivatives), for AggregateBase (when asking for relaxation tensor and RDMatrixPropagator) we achieve better performance for the use. In fact, wanting to work with Redfield tensor in tensor form will be rather a speciality than a norm.

save, savedata, load, loaddata methods for spectra and containers

methods save and load should save and load the whole object with its details (perhaps load should be create an instance of the object and use its _load() method, which fills it? - in both cases an init method to create an empty object is needed). Methods savedata and loaddata, just save and read data from an into the existing object.

Multiple saves of the same object should be a avoided

Saveable still saves the same object multiple times if it is found repeatedly in an object being saved. This has to be avoided.

Create acceptance tests for the mutual support of CorrelationFunction and SpectralDensity

CorrelationFunction and SpectralDensity are closely related. Tests have to make sure that they are mutually compatible.

Create tests for combined Redfield-Foerster calculations

Migrate Makefile to paver

Some of the useful admin functionality which is now in Makefiles should be implemented in pavement.py. In this way the experience a developer has in Windows and MacOS/Linux is the same.

Polynomial approach to the solution of the quantum master equation

Save version of Quantarhei on the top level of the file when saving an object

For reference, a saved object has to know about the version of quantarhei it was saved with

vsignatures #n_particle_approximation

If want to do for example 5 particle approximation for
Vib=[2,3,3]
I think return with no approximation as 5> max(Vin)
But return ndindex for [5,5,5]
I think comparison inside but be (less than )
< Instead of. >

Implement Plottable class

Use it in DFunction, CorrelationFunction, SpectralDensity and AbsSpectrum and other classes to control plotting.

measure_reorganization_energy always returns energy in rad/fs units.

This is related to the units control of correlation function and in general to units control of things that are not directly energy but, say, energy^2

Displacement operator defined differently from standard literature

As pointed out by Ahmed Gharib, the displacement operator returned by the operator_factor class from the package quantarhei.qm.scillators.ho shift by the factor of square root of 2 differently from what is usual in the literature. Full explanation as to why this is so is now part of the documentation, and it can be found here:

https://quantarhei.readthedocs.io/en/latest/classes/qm/ho.html

A simple summary is that the definition of shift operator/displacement operator used in Quantarhei is more convenient in molecular physics. The shift parameter \alpha corresponds to a shift of a coordinate q -> q-\alpha, in a standard definition, one thinks more in a way of "shifting" the creation operator D^{\dagger}{\alpha} a D{\alpha} = a + \alpha.

If no comments appear for 48 hours, I will close this issue.

save and load methods for Aggregate

In this issue, we need to solve the general mechanism how to save objects which contain other objects. It seems that the file-like structure of the hdm5 file should be sufficient for that.

DFunction's at(x) method does not accept array arguments

Let this be considered a bug.