Users.Matlab.Main

About the use of Matlab/Octave in WavePacket

Despite the limited availability of the commercial Matlab software in some academic institutions, we chose that development environment because it offers several useful features. Matlab is rather intuitive to use, easy to learn and offers frequently used numerical tasks (e.g. Fourier transforms, matrix operations, ODE solvers, etc) as built-in functions, allowing for fast code development and high code readability. Moreover, (animated) graphical output is readily available to the user of WavePacket, which is highly relevant for teaching as well as for research projects to develop a more intuitive understanding of quantum dynamics. Finally, the Matlab environment comes with an integrated debugger, profiler, etc.

With the advent of version 7 in 2022, most of the WavePacket codes (except qm_matrix) also run under Octave, though with a few minor restrictions in the graphics/animation output, see also our respective blog entry. However, our codes run substantially slower under Octave than under Matlab.

• Currently the developer is using Matlab version R2022a, running under Windows 11, but WavePacket was reported to run also flawlessly under Linux or MacOS. Reports from users, especially concerning incompatibility problems, on the behavior of the WavePackt software package using other Matlab versions or Octave and/or other operating systems are appreciated.

• No additional Matlab toolboxes are required, the only exception being the 'Control System Toolbox' from which the Lyapunov solver (lyap) is called in qm_balance needed for balanced truncation. If you don't have that toolbox available, there is an alternative Lyapunov solver provided by WavePacket which, however, is substantially slower. Alternatively, you may want to use qm_H2model for dimension reduction.

• If you want to run WavePacket on a computing platform where neither Matlab nor Octave is available, you may want to try out the Matlab compiler.

While setting up and developing WavePacket, we have been especially relying on some features offered by the Matlab/Octave programming environment. Here is a list of our favorites:

• Multidimensional arrays for the points of the DVR or FBR grids are constructed from one-dimensional grids using Matlab’s or Octave's ndgrid function. These grids are then also used to represent wavefunctions, potentials, etc. By virtue of this technique, all of the WavePacket codes are (in principle) ready for an arbitrary number of dimensions.

• When dealing with multiple (coupled) Schrödinger equations, the wavefunctions are stored in a Matlab/Octave cell vector and (diabatic) potential representations are stored in a cell matrix. In both cases each of the cell entries, in turn, is a multi-dimensional array, stemming from the DVR of the position coordinates. Note that Matlab/Octave cell arrays offer additional flexibility due to the fact that individual entries of cell arrays may be objects of different classes. This is conveniently used, e.g., if certain potentials are of different type and/or if certain potential couplings are forbidden due to symmetry restrictions. The same applies to dipole moments and polarizabilities used to describe the coupling of quantum systems to external fields, see also below.

• Since the arrival of version 6, the Matlab/Octave code of WavePacket heavily relies on the concept of object-oriented programming. As is explained here, any WavePacket simulation begins with constructing an object of one of the main classes, either of class wave (wavefunction, represented on a grid) or of class traj (classical density, sampled by a swarm of trajectories) or of class ket_rho (state vectors in energy representation), or any of the sub-classes therof. These objects are then used in the main WavePacket routines, e.~g. function qm_propa for propagation in time which obviously relies on the concept of polymorphism.

• The frequent use of classes and objects offers flexibility in the use of WavePacket for different simulation scenarios. In particular, the different model functions used throughout WavePacket are realized as Matlab/Octave classes, see the following list

folder	model functions
+amo	additional multiplicative operators
+dip	permanent dipole moments (1st order interaction with E-fields)
+efi	external electric fields (1st order interaction with E-fields); typically expressed as a sequence of pulses
+dof	discrete variable representations (FBR) and related finite basis represenations (FBR) for each degree of freedom
+init	initial states/wavefunctions
+kin	non-standard kinetic energy operators, e.g. for Jacobi coordinates or triatomics
+nip	negative imaginary potentials used as absorbing boundary conditions
+pol	polarizabilities (2nd order interaction with E-fields)
+pot	potential energy functions
+sbc	system-bath coupling (LvNE only)
+sht	surface hopping trajectory algorithms
+tmp	numerical propagation schemes for TDSE or CLE or QCLE
+vis	graphical output functions

where the corresponding Matlab/Octave class definition files are organized in package folders, i.e. in subdirectories the name of which begins with a + sign. Note that the implementation of the corresponding objects as Matlab/Octave cell vectors or matrices provides additional flexibility. For example, such a cell vector can comprise objects of different classes thus allowing the use of different DVR schemes for different spatial degrees of freedom

space.dof{1} = dof.fft
space.dof{1}.x_min = -7;
....
space.dof{2} = dof.legendre
space.dof{2}.R_0 = 2;
....

Similarly, cell arrays are used to represent the potential energy matrix, where each entry may be an object of a different class, hence a different model function with different set of parameters.

• Finally, it is mentioned that all the class definitions (already more than 100) used throughout WavePacket are realized as handle classes. In Matlab/Octave, handle class constructors return handle objects, i.e., references to the object created. When passing such objects to functions, Matlab/Octave does not have to make copies of the original objects, and functions that modify handle objects passed as input arguments do not have to return them.