Perils of compiling software under Windows

2024-05-12T10:32:00.846000Z

I have been planning for years to write something in praise of vcpkg as making Windows compilation less unbearable. While building a binary package of
Wavepacket recently, I found, however, that there are also other issues when porting software from Unix, hence a shift in focus.

Introduction

Only after programming C++ under Windows for some time did I fully grasp how much of a hacker's system Unix is. You only need to master a few concepts, such as how the linker works or which package to install for the library headers, then adding external dependencies to your code becomes almost trivial, also for others. And things just work; you might spend ages blissfully unaware of problems such as ABI compatibility, symbol versioning, linker maps or the details of the ELF symbol resolution.

This merry life comes to an end as soon as you decide to compile anything under Windows. Fundamentally, this is caused by different requirements and ecosystems.

Unix distributions thrive on Open Source Software (OSS). At their heart, they are a set of build scripts that download and build thousands of software packages. These packages are easily available through some package manager. Because of that, you can usually assume that every competent user has reasonably up-to-date system package versions and can install all dependencies.

Windows, on the other hand, was meant as an operating system for proprietary software. You have some version of Windows, possibly outdated, on which you install possibly outdated proprietary software that may have been built for an earlier version of Windows. Even worse, the proprietary software may contain
proprietary components that are themselves outdates and have been built for yet another version of Windows.

The net result of many of the subsequent design decisions is that building software on Windows is a moderate pain. This is, by the way, not restricted to OSS, commercial software also has to fight with the limitations. Let us have a look at a few consequences.

Problems (and solutions) when compiling under Windows

1. No package management

Under Windows, the user is assumed to buy software and install it through the software's installer. Package management is not offered by the system itself.

If you want to efficiently manage dependencies, your first step should be the use of a package manager. Windows-native managers such as chocolatey do not cover most development needs. For C/C++ development, you want to use a specific package manager that just builds your direct and indirect dependencies with the latest sources and makes them available to your build. From personal experience, I can recommend vcpkg, but I would assume that Conan is an equally good choice.

The general approach of vcpkg is to offer a CMake (and Meson etc.) integration that mitigates many problems. For example, it sets prefix paths such that the vcpkg are easily found by a package search or copies dlls into your output folder (see the next item).

2. The dynamic loader is comparably dumb

If your code depends on other libraries, you need a dynamic loader to make these libraries available to your code at run time. As a rule of thumb, the Windows dynamic loader loads dlls only from the current directory. It does not understand runpaths, nor does Windows have common directories like /usr/lib where the libraries that you need may have been installed by default (slight simplification, but not relevant for your needs). Of course, user overrides for the dynamic loader like LD_LIBRARY_PATH do not exist, either.

The only general solution for his limitation is to copy all dependencies into the directory of an executable. Even if you build a testrunner, you need to make sure that your library under test as well as all its dependencies are copied into the same directory as the testrunner. Welcome to the state of the art for professional software development under Windows!

A common workaround is to build everything into one common output folder. Then you only need to copy the dependencies once and have them available for all programs. However, this requires of course additional build script infrastructure.

3. The dynamic loader may get in the way

There is another problem with the dynamic loader and plugin architectures. You will usually encounter those if you build a (Python / Java / ...) extension module. Again skipping some complexities, the problem goes like this:

The dynamic loader identifies a dll by its filename.
Dll filenames usually have no version number suffix; symlinks and the like also do not exist.
A program can load every dll only once.

Now imagine Python (which depends on sqlite.dll) loading an extension module that brings an own sqlite.dll. In general, his setup will crash because either Python or the extension module gets the wrong sqlite.dll. The same can even happen if two different extension modules need the same dll.

The solution for this problem is called namespacing. If the dynamic loader is too dumb to distinguish different libraries with the same name, you make the library names unique. Hence the extension module would rename its sqlite.dll into something like sqlite_js9x81nsj.dll.

Note that this situation is usually not a problem under Unix. Ignoring some details, each ELF shared object has an own list of dependencies to search for symbols, so a Python extension module would first look into its own sqlite.so before going up the hierarchy and look into the dependencies of Python.

4. Licenses

It should be clear by now that you cannot reasonably expect even advanced users to build your code under Windows. Instead, you have to build the code yourself, and distribute the resulting binary including all external dependencies.

This opens a new can of worms: Licenses. You might not have noticed under Unix because you only ever distribute your own source code. As soon as you distribute your binary, however, you need to fulfill all license obligations of the external dependencies.

Make sure you supply all used licenses with the binary package. And make sure that you fulfill them. While doing so, you may come across enjoyable special cases; for a harmless example the MSVC C++ standard libraries require you to adhere to the export restrictions of the United States. There may be worse conditions.

5. Miscellaneous other issues

Last, but not least, Windows is simply a different platform, so you need to change the programming style.

When compiling a library, the Unix default is to expose all symbols (functions and globals) in the resulting shared object, while under Windows the default is to not expose any symbols. This difference in behavior makes programs difficult to port. You can change the defaults with flags (for example the CMake target property WINDOWS_EXPORT_ALL_SYMBOLS). However, if the code exports data, such as static variables, this workaround does not cut it, and you will need to expicitly declare the export. as I had to do with my underlying tensor library
By default, what is a shared object file "mylib.so" under Unix is split into three different files under Windows: the runtime dependency "bin/mylib.dll", the debugging information "bin/mylib.pdb" and the import library with information for the linker "lib/mylib.lib" For additional confusion, the latter has the same extension ".lib" as a static library.
If you use the Windows API directly, always convert your strings into UTF16 wide strings and use the "W" version of the Windows API functions. This guarantees consistent behavior, the narrow string "Ansi" functions may depend on the user localization. This is a leftover artefact of Windows' early support for Unicode.

Some other issues that I have never come across personally, but that you should be aware of:

Never pass memory ownership over Dll boundaries, allocation and deallocation must always be done in the same module. If you build a library that hands out allocated memory, offer an API to free this memory again. The only exception is if you can guarantee that all affected modules are build against the same runtime. Failure to do so can lead to mystery crashes.
The background is that different modules (possibly clients of your library) can be built against different runtimes that have different ideas about how to allocate memory. This problem is rather esoteric under Unix; by default all shared objects use the same global libc. These differences can be a major pain when porting a library to Windows.
You are more likely to encounter ABI compatibility problems under Windows, because the system has not been compiled with essentially the same compiler. Plus, there are different calling conventions in use. Most likely, you will be fine if you at consume or expose only C interfaces, but be aware
of this issue.

Conclusions

This list of problems is by no means exhaustive, but should cover the main pain points. I hope it gives a reasonable overview of the main pitfalls and enough
pointers and keywords to solve the issue.

Personally, all the problems that I have also seen at work have made me into a very practical advocate of OSS. Once they are packaged, OSS components are easy to use and upgrade, usually make sense as a package, and do not require additional build infrastructure. This is notably different from issues that we have at work with commercial software.

Wavepacket 0.3.6 released

2024-05-09T14:39:21.283000Z

The main focus of this release was an improved compilation. In particular, being able to provide Windows binaries, given that this is a recurring issue.

1. Windows binaries and simpler compilation

Version 0.3.6 is the first to ship a Windows distributable, consisting of a Python package with associated documentation and Demos. Everything is still somewhat crude (documentation, setup, required preparation), but it works at least.

To get there, a host of other changes needed to be done. Wavepacket now depends on an upgraded version of the tensor library (which builds rather easily on Windows). Another side effect was that I had to revisit how to build Wavepacket easily under any operating system. The upgraded Vagrantfile now provides a template for compiling all dependencies that should be rather portable to arbitrary systems.

2. 2D contour plots

To slowly improve the plotting capabilities, this release adds another plotting type. The plot is still somewhat simple compared to the Matlab version, but it is a first step towards more complex plotting setups.

3. A solver based on Faber polynomials

This one has been on my personal agenda for years. I wanted to implement a polynomial solver similar to the Chebychev solver, but suitable for open quantum systems. The new Faber solve is moredifficult to use than the Chebychev solver because we need to specify a domain of convergence in the complex plane, not just an interval on the real/imaginary axis. However, the final result should be reasonably straight forward to use. See https://wavepacket.sourceforge.io/cpp-doc/current/demo_faber.html for the ugly details.

Because most of the effort in using the Faber solver is spent on figuring out the support of the Hamiltonian under study, this issue also spawned some additional functionality to get extremal potential values, to truncate simple operators (DvrOperator / FbrOperator, which are diagonal in the DVR or FBR), and to estimate the maximum eigenvalue with the power method.

4. Smaller features

As usual, other issues got fixed along the way I would just refer to the detailed changelog under https://sourceforge.net/p/wavepacket/cpp/git/ci/master/tree/NEWS instead of repeating it here.

Outlook

I currently plan to finish the API, which will introduce some breaking changes over the next releases, and to implement the final functionality to be roughly en par with the Matlab version in terms of basic functionality. I also have a mild interest in fleshing out the Python interface more, because Python offers a more natural trial&error approach to running Wavepacket simulations.

Wavepacket 0.3.5 released

2023-04-10T21:05:58.021000Z

This release took unusually long; I usually strive for about two releases per year. This one was overdue due to some assistance with getting the underlying tensor library to use CMake, and a rather discouraging set of issues that I mostly dropped for an alternative approach. Anyway. What did change?

1. Much better support for open-system dynamics

Wavepacket can propagate a density matrix with some system-bath coupling completely in the DVR. That is pretty neat, but, as it turned out, mostly useless.

When you use a reduced density matrix for open system dynamics, you always set up your Lindblad or Redfield Liouvillian in the energy eigenstate basis. This has two reasons. One is that you generally go for the Caldeira-Leggett model where the bath consists of uncoupled harmonic oscillators at different frequencies, which defines the coupling in an energy eigenstate basis. The other reason is that you normally want to have relaxation to thermal equlibrium as limiting case, which is difficult to express outside of an energy eigenstate basis.

Version 0.3.5 offers functionality where you only need a system Hamiltonian and specify a spectral density and get a Lindblad Liouvillian with a moderate amount of work. See the example MolVibration/OH/3 for a slightly more involved example of that. This scheme, however, does not deal with the Lamb shift, which can be relevant, but is tricky to calculate in any case, so I will probably stop at the current level of assistance.

2. More consistent names

I was not happy with some of the names, sometimes even order of arguments. So there was a first round of API-breaking changes: consistent naming and ordering.

For the complete list of breaking changes, see the NEWS file. Just to highlight a few:

"PropagatorPrimitives" are now called "Solvers", so for example ChebychevPrimitive => ChebychevSolver
dropped useless prefixes and suffixes. For example, "Interval1D" is now called "Interval", or "PlotObserver" in the Python module "wavepacket.plot" is now called "Observer" (the module name says all about plotting)
switched the order of arguments in the constructor of "State": first grid, then the coefficients. That is the order used everywhere else, and I always mistyped that.

3. Build cleanups

I went over the (CMake) build and cleaned it up in a lot of places. There should be fewer test errors, generally more stability, a few less dependencies etc.

The most noticeable change is that I now try to use the installed googletest version before downloading something.

4. Smaller features

As usual, some things cropped up along the way:

new functors: Lorentzian and SoftRectangularFunction
function to get the lowest N functions from an eigenstate solver, function to orthonormalize a set of states
function to construct the adjoint of an operator
made Python demos always produce plot output where possible

Outook

At work, one of highlights in the last year has been to add a Python interface to some code. As a consequence, I am now way more familiar with Python than before, and would like to redo the Python integration. Also, the underlying tensor library has been adapted to CMake, which makes it way easier to compile, but requires code adjustments. Put together, I hope to finally provide a Python wheel, also for Windows, for the next release.

Breaking changes will continue, there is one large set of changes (Python is very opinionated about snake_case), and a few smaller issues with consistency, but I hope to get over these within the next two releases and finally keep the API mostly stable.

Otherwise, I am going to add a few more usability goodies. Wavepacket is already pretty usable, except for being slow when doing open systems and poor in the plotting department, so these would be the other focus areas.

Wavepacket 0.3.4 released

2022-01-11T20:13:23.124000Z

New features in this version:

There is now a much simpler approach to build Wavepacket with a
virtual machine.

Install two programs (Vagrant + Virtualbox), modify a config file,
open a console, type three to four commands, and you can execute
Wavepacket Python code in your Web browser. For more details, see
README_VirtualMachine in the root directory. Needs some more
polishing, but should be reasonably plug&play.
Plotting in Python was improved considerably.

Plots are still limited to one dimension right now, but you can set
the units to use for coordinates, energies etc., and it is possible
to save the plots as animations.
Added exponentiation of operators.

This feature is still in its infancy, but at least basic
functionality for certain operators is there. This will allow
absorbing boundaries for the efficient Chebychev propagators or the
split operator method.

I also removed some temporaries to speed up the Chebychev
propagators, and did some cleanup that will help with a Windows
release.

Dependency management with DotNet or Why you should use Paket

2021-09-09T21:17:20.795000Z

Today something completely different. I wanted to write about this topic for a long time, and it also has nice overlap with vcpkg (which I also wanted to write up for some time now), so I finally decided to just do it.

This blog post was inspired by me learning more DotNet internals than I would normally care about. One of the rather frustrating experiences during this process was the lack of readily available information on how DotNet deals with dependencies. Hence this blog post.

On second thought, maybe the information was there, but perfectionist that I am it took me some time to accept the state of the art.

One word of warning: I have setup our DotNet code at work to use Paket for dependency management and never again looked back to Nuget. It might be that the situation discussed here has improved by now.

The whole article is licensed under the CC0 license. Copy parts, whatever. It would be nice to receive attribution, but in the end I do not care.

Shortest possible summary

The three takeaway messages (and sections) of this post are:

Dependency management is hard.
Be careful and make conscious decisions, or you will be hurt in places you did not even know.
Dependency management in DotNet is harder.
You are limited in your options, and it is rather easy to break things.
Ditch Nuget (the tool), choose Paket instead.
Nuget is complicated, lacks functionality, and is terribly brittle in edge cases.

1. Dependency management is hard

When you have a project, you typically decompose it into coherent modules. These can be executables, test runners, libraries etc. The modules expose a public API that allows other, consuming modules to use the encapsulated functionality. This public API changes over time, either in syntax or in behavior. Hence, modules always reference other modules by a specific version, either implicitly (e.g. "the version from this commit") or explicitly through some version number, e.g., from semantic versioning.

In almost all non-trivial cases, not all code will be written by you / your team, and you will rely on external modules for certain functionality. In this case, you need to make sure that these provide the correct API version for error-free consumption by your code. That is, you need to manage external dependencies of your code.

Note that in the following, I will switch to VisualStudio / DotNet speak. VisualStudio calls projects "solutions", and the individual modules "projects" with "assemblies" as build output. Finally, assemblies can be wrapped in a (usually Nuget) "package" that contains additional metadata, in particular dependencies on other packages.

Example 1

Consider the trivial case of a program P depending on an external assembly E, as depicted in figure 1

Figure 1: simple dependency of a program P on an external assembly E in a specific version.

All you need to do here is to choose the correct version of E and you are done.

The dependency is evaluated in two contexts: During the build, and at run time. During the build, the compiler looks for the assembly E and essentially verifies that all referenced symbols with the correct signature exist. When running the program, the DotNet runtime loads the assembly E with the correct version when needed and redirects calls to the corresponding symbols in E. In-between, the symbols are stored essentially as strings.

Example 2

Let us choose a more complex example, where P depends on two assemblies A, B, and these in turn depend on an external assembly E. See figure 2 for the scheme.

Figure2: Fundamental diamond dependency conflict, where the program P depends on two assemblies A, B that reference an external dependency in different versions x.y.z and a.b.c, respectively.

In practical situations, the dependency graph can be more complex, with transitive assemblies in between or more than two assemblies requiring E, but the circular structure as central feature is commonly found.

To understand the principal problem, consider the following scenario. Let A be the output of an internal project that is under your control, while B is an external assembly that you only get in compiled form. You compile A against E in version x.y.z, and at runtime you also only provide E in version x.y.z. to fulfill both dependencies from A and B. Let the versions x.y.z and a.b.c be different.

Now assembly B calls a particular API (E in version a.b.c) but gets a different API (E in version x.y.z). In particular, symbols or the underlying behavior may have changed between the two versions. Whenever B calls one of the changed symbols in E, bad things may happen. In the best case function signatures changed, then the runtime cannot find the requested symbol and throws an exception. This can cripple your application because some functionality is not reachable or even crashes the application, but at least you can localize the error. In the worse case E behaves differently from what B expects. Then you have a very obscure, hard to trace down bug.

What is particularly insidious in this case is that all the code is perfectly sound. No amount of static code analysis and almost no review will highlight an error. However, your total system, composed of P, A, B, E(x.y.z) is broken, and this only at runtime. In the literature you will find such situations labelled with the term "hell", as in "dll hell" or "jar hell" or, as I will expand below, "assembly hell".

What to do with a diamond dependency

There are basically three approaches:

Pick & Pray
The simplest solution uses a common version to satisfy both dependencies. However, as discussed above, this is largely programming by coincidence and not guaranteed to work.

There is one exception, though. If and only if the assembly E is trustworthy not to change the public API without notice, and if no such notice has been given between x.y.z and a.b.c, the two versions should be interchangeable. For example, if E uses semantic versioning and the two versions differ only in the minor or micro version number, then you can pick the higher version to fulfill both dependencies.
Duplicate the assembly

Another way to solve the diamond dependency is to resolve the diamond like in figure 3.

Figure 3: Possible solution to the dependency problem of figure 2, by supplying assembly E twice with two distinct versions.

If you depend on the assembly E in two different versions, well, you just supply the assembly in two different versions. This generally works, however, there are a few caveats.

If the assembly E has internal state, this state will be duplicated, which can be a problem. For example if E is a logging framework, duplicating E gives you two logging frameworks that work independently. These can then interfere in arbitrary ways, for example by overwriting each other's files.

More importantly, this solution is only available if the runtime supports it. Some runtimes that do not support this setup are the ELF format (linux shared libraries) or Sun Java before version 9 (though you can work around both restrictions with modest extra effort). In both cases, the application holds a list of modules, either jar files or shared libraries, that can be searched for missing symbols (functions or variables). When a symbol is requested, the runtime or dynamic loader takes the symbol name essentially as a string, and goes through the list of modules until it finds a matching entry. If multiple entries exist, only the first is ever returned, hence duplication does not work.
Enforce a common assembly version

You can try to ensure a.b.c == x.y.z, that is, that all code only ever requests a single version of the external assembly. This should be your preferred solution whenever possible.

Of course, as usual, this approach is not always feasible. It only works if you have complete control over the affected assemblies A, B in figure 2. But often, at least one of the assemblies is itself an external assembly; then you would need to modify and compile foreign code, which is considerable extra work and complexity. Even if the assemblies A, B are in principle under your control, they might be maintained by different teams. Or they may be used by several internal clients with different requirements. Requiring a consistent version of all dependencies can then add a lot of inertia to the development process. Instead of a "hey, let's quickly update this package" approach, you need to involve and synchronize multiple teams etc., which can also translate into considerable costs and effort.

2. Dependency management in DotNet is harder

Now how does the DotNet environment fare with regard to the dependency problem? Unfortunately not too well.

Apparently the runtime seems to support the layout of figure 3. You can have the same symbol name referring to different code in different assemblies. However, the assembly loader effectively does not support this layout. The details are listed here, but in essence every assembly has a name (and version, but this can be overridden), and assemblies are referred to by their name. The DotNet assembly loader can load an assembly with a given name only once, unless you use rather arcane features like loader contexts.

In essence, unless you go quite some extra miles with renaming assemblies, rebuilding other assemblies to refer to the new names etc., the duplication of assemblies does not work. You may spend the extra effort for one or two crucial assemblies that have been shown to be problematic, but not for the dozens of dependencies that a serious DotNet application may depend on. So in general you can choose: Try to get common dependency versions or pray.

While DotNet is not alone with this problem, its impact is exacerbated by the fact that dependency trees can quickly become deep and wide. In particular Nuget packages like the various Azure components have a tendency to drag in the whole world in arbitrary versions.

To give some numbers from work: We have a simulation, probably 50-100k lines for the kernel, that allows simulation (and corresponding tracking for billing, authentication etc.) in Azure, which is why it depends on various Azure packages. The program has more than 200 external dependencies. In contrast a multi-million-line, highly complex desktop simulation has less than 50 dependencies, many of them highly specialized and conflict-free.

For the versioning problem, a notable example is Newtonsoft.Json, probably the most widely used DotNet Json parser.
Various Azure packages depend on Newtonsoft.Json in version 6.x, while other packages required it in a version like 12.x. That is a whooping difference of six in the major version! It seems safe to assume a priori that these dependencies are incompatible.

A third problem comes from the way Nuget packages tend to deal with version constraints. The default version constraint, and what feels like the only one I have ever seen is a minimum version. So if you have a Nuget package/assembly X that has been built against a package/assembly Y in version, say, 2.1.0, then the Nuget package for X declares a dependency on the package Y with any version >= 2.1.0. This dependency can be fulfilled by 2.1.1 (ok), 2.5.0 (probably ok), but also 10.0 (not ok). While this approach has its virtues, it makes dependency management harder because you lack warnings that would indicate potentially incompatible versions.

3. Ditch Nuget (the tool), choose Paket instead.

Given the state of the art of DotNet dependency management, there are a few features that a DotNet package manager should have:

When you have a solution, consisting of libraries, test runners, programs etc., you generally want a single version of each external assembly for everything.

This is pretty important. Consider a test runner using a different version of some external assembly from your production program. You may have, for example, a dependency error in your program that is never caught by a test or different behavior in test and production. Good luck debugging that.
You should be able to get a dependency graph for your solution.

The graph shows all dependencies from your projects to Nuget packages and transitive dependencies between Nuget packages. It allows you to review the versions of all dependencies, so that you can identify potentially problematic dependencies. Once you have "signed off" the dependency graph, it should be the immutable reference for the packages/assemblies to use, and should only change if you add/remove or update dependencies.

Unfortunately, Nuget has the very fundamental design decision; Nuget manages dependencies project-wise, not solution-wise. In theory, this offers additional flexibility, and is probably the only way to integrate with the build system. However, it adds a lot of overhead to the dependency management, because you need to manually ensure that all projects use the same version of a given package. I have fond memories of a coworker spending hours and an endless stream of curses on this task for a comparably trivial program (~50k lines), and people fearing an upgrade because of the days of trouble it brings. Also, while Nuget nowadays offers something like a dependency graph, you get one for each project, which limits its use.

Nuget has a few further downsides:

No C++/CLI support.

If you combine DotNet and native code, you might sometimes like to consume Nuget packages in the C++/CLI glue code. Unfortunately, Nuget has no C++/CLI support whatsoever. C++/CLI projects cannot have Nuget-managed dependencies, and when looking for potential incompatibilities between packages, Nuget never looks past a C++/CLI project. That is, to spell it out, you do not even get a warning if your dependencies are entirely broken.
Integration into the build system.

By default, when you use Nuget, it integrates into MSBuild. Superficially, that sounds pretty cool, and allows, for example, to just specify the Nuget packages in your project files. In practice, it is not.

MSBuild is simply a pretty complex monster, and bolting down yet another layer of functionality does not make it more robust. I would recommend to turn on detailed logging for a project with some dependencies, and read up where and how MSBuild looks for potential assemblies to link your code against.The result should convince you that you want to override MSBuild, not integrate. There are also some exotic cases, where MSBuild can "lose" assemblies in interesting ways that suggest you do not want to integrate.

Entering Paket, an alternative package manager. It is not perfect: the documentation took some time to get used to, and the way it manipulates project files feels a bit dirty. But Paket supports the work flow as it should be, and it has this down-to-earth feeling of Unix tools: doing one thing properly.

To use Paket, you logically proceed in two steps. First, you define what external dependencies with which version constraints you need globally for your solution, and run Paket on this input. The result is a dependency graph in a file called paket.lock that tells you every single package and version that is needed. This in itself is already pretty cool: Just by parsing this file, you can see which package requires what other package. If you are not happy with the dependency graph, you can fine-tune by adding more requirements or changing required versions.

Once Paket has created the dependency graph it avoids touching it. When you add or remove dependencies, it will build a new dependency graph using the existing one as starting point, and you can also update packages, which apparently changes dependencies. But as a general rule, once you have created a dependency graph and are happy with that, you never need to worry about breaking something by accident. Afterwards, you need to add a pre-build step that downloads all the packages in the dependency graph and puts them in some specific directory.

The second logical step is defining for each project what dependencies are needed. Paket will take a look at these definitions and modify the VisualStudio project files. For every dependency, it adds an <assemblyreference> tag to the project file with the location of the downloaded package as HintPath. This is about the strongest coupling you can have to an external assembly in MSBuild. MSBuild will always look in the HintPath first, and copy the assembly to the output directory in any case etc.</assemblyreference>

That is all that Paket does. No complicated project manipulations, no complex interaction with the build system, just downloading packages and forcing dependencies to the downloads. And assembly references also work perfectly well with C++/CLI, so with Paket you can add Nuget packages to C++/CLI projects.

Wavepacket 0.3.3 released

2021-08-08T11:14:46.870000Z

The version features two major improvements towards the goal of being more usable:

The Python interface has a plotting observer.
There is now an Observer (written in Python and used in the C++ code) that utilizes Python's matplotlib to plot the wavefunction on the screen. This is currently not that impressive because it works only for one-dimensional wave functions . However, with the technological barriers now overcome, you can expect more plotting options for the Python interface in the future.
See the DoubleWell Python demo for an example.
The operator hierarchy has been rectified a little.
This may give a minor speed boost in certain situations, because real-valued and complex-valued operators are now disentangled and because real algebra is faster than complex algebra. The main gain is that it should be easier to reason about most operators in code, which is a prerequisite for things like better optimization of operators, or the split-operator scheme

Besides these two major improvements there have been a couple of smaller usability or special-usage improvements. The most noticable:

Removed the OdeRelax propagator, since it should be always slower than the Chebychev relaxation; the latter can also propagate in negative imaginary time now.
Diagonalization always returns explicitly real-values results.
You no longer need to construct a projection operator to project() onto a coupled channel; you can now simply give the index of the coupled channel.
Most Python functions no longer need casting between Numpy arrays and the C++ tensors.
ProjectionND can take a state as input to project onto, not a tensor + grid (which defines a State)

Unfortunately, I failed another time at providing a Windows package, but there will be a much more convenient Virtual machine with the next release.

Convergence #2: time steps

2021-04-19T21:15:32.585000Z

Note: This post should probably be merged with the Wiki, in particular the page on the time-dependent Schödinger equation. The second part should be easy, but the first part uses a pretty different angle, so I first post the whole thing here.

This document describes how to propagate a wavefunction in time. The question should be answered are:

What time step should you use?
How can you monitor during a propagation that your chosen time step is ok?
What are advantages and disadvantages of different propagation schemes and when should you choose which?

To restrict the scope, I mostly skip convergence issues with respect to the grid. So either you already converged the grid, or I just take Hamiltonian plus numerical errors from the finite grid as The System and disregard deviations from reality until they come back for revenge.

Furthermore, I will only consider time-independent systems here. Basically, if you can converge a time-independent system, you can also converge a time-dependent system: Guess a slightly smaller time step and adjust until you are done. However, introducing the required formalities for proper reasoning (cavity-dressed states etc.) would bloat this already long description even more.

Some theory

Model system

As a model system throughout the text, I take one of the Wavepacket demos as blueprint. The demo, HarmOscillator/Gaussian_1D/1, describes a one-dimensional harmonic oscillator (omega = m = 1 a.u.), where the initial state is a displaced Gaussian wavepacket.
The native time scale is the classical oscillation period, which is a.u. for the given parameters. Apart from an offset of 1/2, the eigenenergies are spaced by the energy omega, hence the time scale for the n-th quantum eigenstates is approximately the n-th fraction of the classical oscillation period.

I keep the total simulation time to about one classical period. Longer times would make the simulations slightly slower, but do not contribute qualitatively new features.

Figure 1: Projection of the initial state onto the energy eigenstates as a function of the eigenstate's energy. Red crosses: 128 grid points, interval [-10,10] a.u.. Blue dots: 512 points, interval [-20,20] a.u.

The first figure shows the projection of the initial state onto the eigenstates of the harmonic oscillators, that is, their populations as red crosses. The Python script that generated this plot can be found here. The populations decay roughly exponentially down to about 10^-9 at an energy of 50 a.u. and remain on this level. This latter feature is not what we would expect from an analytic solution, instead the populations should decrease monotonically with energy. In fact, this feature is an artefact from the finite grid, which is particularly pronounced for high energies. To demonstrate this, we can use a brute force approach and double the size and density of the grid (increasing the number of points fourfold). In this case (blue dotted line), the populations decrease further as expected.

This observation highlights two issues that we will encounter again further below:

We have an unnaturally high population of high-energy states. This will cause trouble during propagation.
We can safely truncate the energies of the Hamiltonian at 50 a.u. This will perturb all higher energy states, but they are garbage anyway.

As a final remark, this system lacks some of the real-world complications for simplicity. The Hamiltonian is Hermitian and time-independent, there are neither open-system terms and no absorbing boundary conditions. We will briefly touch some of these complications later when discussing different solution schemes.

Finding a good time step and monitoring convergence

A first intuitive approach to determine a good time step is based on a qualitative understanding of the system dynamics. If you plot the spectrum of your initial state, as in figure 1, the fastest dynamics can be expected to come from populated eigenstates with the highest energies. So if you choose a time step that is considerably smaller than the oscillation, say one tenth, you would expect the simulation to converge. For the harmonic oscillator example, you could generously say that states above 50 a.u. have negligible population.
An energy of 50 a.u. corresponds to an oscillation time of . Taking one tenth of that, you could reasonably guess that a timestep of 0.01 a.u. should be good enough for reasonably accurate numerical schemes.

This approach is certainly useful to get a good order of magnitude guess, however, it is also not correct. We will go into the details when we examine our example system in the next section.

Another standard approach is to try out different parameters, that is, time steps, and to check how much the results differ as a function of the time steps. The differences should become smaller as we go to smaller time steps, i.e., more accurate solutions. Hence, we only need to set some qualitative target when we are satisfied, and run a couple of calculations.

As a simple measure for convergence, we can take the induced metric of the L2 scalar product. If v_1 and v_2 are two results with different time steps, then their distance is

Here, I assume that the norm of the states is approximately one. Since only relative quantities and order-of-magnitude estimates are important here, I also dropped the prefactor sqrt{2} for simplicity.

However, there is also a much simpler metric available. For this, we first note that in general the norm of a wavefunction is constant over time. This is a direct consequence of the time evolution operator being unitary. All numerical schemes formally produce approximations to the time evolution operator, and unless they are carefully crafted, these approximations are not unitary. Hence, most numerical schemes fail to preserve the norm, and we can easily find an ok time step by checking that the norm of our wavepacket does not change.

However, like all good ideas, there are certain cases where this scheme fails. In particular, as soon as we use absorbing boundary conditions, the norm ceases to be a preserved quantity.

Application to model system

Now let us see what happens when we play around with the harmonic oscillator model system. Figure 2 shows the norm of the wavefunction as a function of time for different time steps ( The Python script can be found here). We used essentially a fifth-order Runge-Kutta solver to evolve the wave function in time.

Figure 2: *Norm of the wavepacket as a function of time for different time steps.

Recall that we estimated a good time step as 0.01 a.u. For a too large time step of 0.1 a.u., the norm diverges essentially immediately, showing the immediate use of the norm as a monitoring tool. However, even for our estimate of 0.01 a.u., the norm diverges rapidly at some point, while a slightly smaller time step of 0.005 a.u. shows no divergence.

So, what has happened here?

dt = 0.01 a.u.	dt = 0.005 a.u.

Figure 3: *Projection of the time-dependent state onto eigenstates of the Hamiltonian, similar to figure 1, for different times. Left plot: dt = 0.01. Right plot: dt = 0.005.

To gain some further insight, we calculate spectra similar to figure 1, but for different times. This is done in figure 3 for the two relevant time steps of 0.01 a.u and 0.005 a.u. (Python script here).

These spectra tell a story. Apparently, the numerical errors that lead to the norm deviation accumulate not in the energy range where most of the wavepacket resides (< 50 a.u.), but at the highest energies. Now if we have an eigenstate for a given energy, we try to approximate the correct time evolution, exp(-iEt), by discrete time steps, which gives an error. An example is shown below for the exceedingly bad but simple Euler method. Fifth-order Runge-Kutta methods are qualitatively similar, although their error is much smaller. In particular, you need to sample one oscillation by multiple steps, say 5 or so, otherwise, you cannot hope to describe the dynamics accurately.

Figure 3b: *Approximation of the correct time evolution (black circle) with the Euler method for different time steps as fractions of the oscillation time T.

Now for the highest energy states, we have E about 200 a.u., and hence an oscillation period of about 0.03 a.u. Such an oscillation is poorly sampled by the time step of 0.01 a.u., leading to a rather large error in the norm. For the half as large time step, we sample an oscillation with about 6 discrete time steps, which is just enough for reasonably accurate results. Note that if you look closely, the highest energies also show a norm increase for dt=0.005 a.u., but much smaller, maybe one order of magnitude every 10 a.u. Hence, if we would simulate much longer, say, for 100 a.u., the norm would also diverge for the smaller timestep.

Corollary

The main result from trying out the theory is something like

Not the energy / dynamics of your initial state dictate the required time step, but the whole spectrum of your system.

This observation has an important consequence: If you are interested in performance, that is, large time steps, you should always truncate the spectrum of your Hamiltonian or its components! There are two reasons why truncating operators should not bother you too much:

On the purely numerical side, these high-energy states are often garbage anyway, as discussed in figure 1. Hence, there is little harm done if you screw them up completely by truncating the operators.
On the physics side, you need to recall that most systems are approximations. For our harmonic oscillator example, we observed severe numerical errors at energies larger than E = 50 a.u., which corresponds to the fiftieth harmonic oscillator eigenstate. Recall that harmonic oscillators are low-order Taylor expansions of the true potential. So whatever the real system is, we can safely assume that, say, the fiftieth excited state is not going to look like a harmonic oscillator eigenstate.

No matter how accurately you model your system, you will always have such severe approximations. Maybe you use the Born-Oppenheimer approximation and ignore non-adiabatic couplings: This approximation fails for large energies. Maybe you use a model for your thermal bath: Rest assured that your model is not quantitatively correct at high energies.

So: always know and pick the energy range that you consider important. Then be brave and truncate the spectrum of your operators (typically kinetic energy and potential) outside of this range, knowing full well that you limit the validity of your calculations to sufficiently low energies. What you loose in correctness, you will definitely gain in computation speed and numerical robustness.

Figure 4: *Norm as a function of time for different time steps. The kinetic and potential energy have been truncated at 50 a.u. each, that is, values larger than 50 a.u. have been set to exactly 50 a.u. in the corresponding DVR/FBR grid representations.

For the model system, figure 4 shows the effect of truncating the kinetic energy and potential to 50 a.u. each (Python script here).
We can clearly see that the norm no longer diverges for the time step of 0.01 a.u., giving us a factor of 2 decrease in computation cost.
In practice we might truncate at considerably lower energies (20 a.u. or so); here we chose 50 a.u. because that is where the obvious numerical artefacts start.

Solution schemes

To solve the time-dependent Schroedinger equation or Liouville-von-Neumann equation, different methods or method classes exist. In the following we will present some important solvers and discuss their advantages and disadvantages

Runge-Kutta and other ODE solvers

The basic idea for a Runge-Kutta or similar scheme is to expand the time-dependent wavefunction in some orthonormal basis

Note that the grid representation is such an expansion with the coefficients c_n(t) being the values of the wavefunction at the grid points psi(x_n, t). If we plug this expansion into the time-dependent Schroedinger equation, we get

The latter is a set of ordinary differential equations (ODEs) coupled by the matrix elements of the Hamiltonian.

There is a long history of numerical approaches to solve ODEs, and many solvers exist. Prominent general-purpose solvers are the Runge-Kutta schemes, the Bulirsch-Stoer method or various multistep methods.

As a particular case, there are also adaptive solvers, for example adaptive Runge-Kutta methods. Here, a Runge-Kutta method of order N gives you the solution for a given timestep, then you check with a special other method of order N+1 if the solution looks reasonably converged. If the errors are outside of a given bound, the solution is recalculated with a smaller time step. If, on the other hand, the solution seems extremely well converged, the time step can be increased instead.

Typical parameters that you supply are the order N of the scheme, and the timestep \Delta t. Alternatively, adaptive methods require an error bound, usually consisting of a relative and an absolute error. Generally, for a (Runge-Kutta) scheme of order N, the error scales as .

Advantages

ODE solvers are incredibly robust work horses. They work with time-independent problems, time-dependent problems, complex Hamiltonians (absorbing boundary conditions), open systems, ... Just throw a problem at them, and you get a reasonable solution.

If you take an adaptive solver, you do not even need to bother with the timestep directly. Just supply an error bound, and the adaptive solver spits out a reasonable result.

In short: ODE solvers tend to always work, and they are very easy to use.

Disadvantages

ODE solvers do not use any knowledge about your system or quantum mechanics in general. This is boon because they are robust, but also a bane because it makes them the slowest of all solvers. There are also some cases where convergence is poor, which again translates into small time steps and low performance.

In short: If computation time is an issue, you should generally avoid ODE solvers.

Split-operator method

The basic idea of the Split operator method is to separate the Hamiltonian into components that you can diagonalize easily,

If you can diagonalize an operator, that is decompose it into a basis transformation T and a diagonal operator D in this basis, , you can also easily exponentiate it as

as can be seen by writing out the Taylor series of the exponential function.

The split operator then consists of writing the short-term propagator as a product of individual terms. For the common Strang splitting, we get

Note that the right-hand side is by construction unitary. As a consequence, the split operator method preserves the norm of the wavefunction.

Besides the time step, there are two further "parameters" to choose. One is the order, though in practice, everything beyond the Strang splitting tends to become complex. The other is the choice of operators H_i; again, in practice, this is usually intuitive. You normally choose the potential as one component (diagonal in the grid representation / DVR), the kinetic energy as the second (diagonal in the underlying basis / FBR), and some dipole couplings or so as the third one for time-dependent problems.

Advantages

The split-operator method has a number of unique advantages.

First of all, despite the rather low order of , it is generally rather accurate. As a handwaving argument, you can consider the artificial case where only the Hamiltonian is chosen and diagonalized. Then the split-operator method is actually the exact time propagation. As an alternative hand-waving argument, the conservation of the norm means that the approximate solution is searched only in a reduced-dimensional manifold of states with unit norm, which gives more accurate results.

A second advantage is that the split-operator method with Strang splitting is rather cheap for simple systems. Consider a Hamiltonian decomposable into two terms, kinetic and potential energy. Consider further that you are generally not interested in the value of the wavefunction at the end of each small timestep, but only every hundredth or so timestep. Then you can combine the last term of the n-th timestep and the first term of the (n+1)-th timestep into one exponential term. As the exponentiation can be performed beforehand, you need to apply two diagonal operators and perform two FBR <-> DVR transformations for propagating a single timestep \Delta t. Contrast this with a second-order Runge-Kutta-scheme with similar error of . There you would have to apply the Hamiltonian twice for a total of four diagonal operators and four FBR <-> DVR transformations.

The split-operator scheme can also be extended to time-dependent systems, but becomes less efficient there. In practice, you then have to exponentiate a diagonal operator for every small time step.

Disadvantages

The split-operator method cannot be applied if your system has terms that cannot be easily diagonalized. This prevents its use with certain Hamiltonians, for example products of momentum and coordinate operators. Also, open-system terms tend not to have an accessible diagonal representation. Hence, although you could in principle generalize the split-operator method to density operators (replace Hamiltonians by Liouvillians), this is not very useful in practice.

Another issue concerns the implementation. An efficient diagonalization of certain operators like dipole couplings is possible, but requires some fiddling. This is, for example, the reason that the method is not yet implemented in the C++ version of Wavepacket as of version 0.3.2.

Polynomial methods

The basic idea of the various polynomial methods is an expansion of the short-time propagator into some polynomials. In the simplest form, this can be written as

although in some schemes the coefficients also depend on \Delta t. In particular, for various classical polynomials, the P_n can be calculated through a recursion relation, and the coefficients c_n can be precomputed. In practice, you then apply the Hamiltonian N times with some simple algebra in-between.

The free parameters are the choice of polynomials, the order of the expansion N, and the time step. The error of an N-th order expansion is generally . Polynomials are typically chosen based on some additional properties.

Particularly useful are Chebychev polynomials. These have the property that they converge continuously against the exponential function within a certain domain. The details are a bit involved, for example you have to normalize the spectrum of your Hamiltonian to the range [-1, 1]. However, as this is implemented in Wavepacket, you just need to specify a maximum error, from which the order of the expansion is calculated automatically. The time step for the Chebychev expansion is inversely proportional to the spectrum of your Hamiltonian, hence this method benefits a lot from aggressive truncation.

Note that the order of expansion should generally be rather large (N > 20) for this method to be efficient. However, if you go to too large values (N > 100), you can experience numerical instabilities. Basically, the polynomials themselves have huge function values that can even overflow the double-precision range. So you should try to stay with not too large polynomials, say N < 60.

For these large orders, typically, you will choose a timestep that is much larger than for, say, an ODE solver. Often, this is not a problem, as the numerical time steps tend to be much smaller than the usual dynamics. However, if you do need the wavefunction at high time resolution, you can generally get there with additional interpolation steps. These require some work to set up, but are very efficient; in particular you should not need to apply the Hamiltonian to a wavefunction for the interpolation to work.

Advantages

The huge advantage of polynomial expansions is their efficiency.

To get an estimate, let us consider the Runge-Kutta procedure with the truncated Hamiltonian of figure 4. At a timestep of 0.01 a.u., and with the fifth-order Runge-Kutta procedure, we need to apply the Hamiltonian 5 x 100 times to propagate the wavefunction for 1 a.u., with a corresponding number of vector additions, multiplications etc. In contrast, if we take an expansion in Chebychev polynomials with a spectral range of 50 a.u., a timestep of 1 a.u., and an accuracy of 10^-8 ( about the order of the numerical errors, see figure 1), we need just 46 polynomials. That is, to propagate for 1 a.u., we need to apply the Hamiltonian about 50 times and require only one tenth of the computation time for the Runge-Kutta method.

Disadvantages

Polynomial methods in general, and the Chebychev expansion in particular, also have severe downsides.

Most importantly, all polynomial methods fail for time-dependent Hamiltonians. The reason is that the simple Taylor series for the short-time propagator turns into a hideous expression with time-integrations of time-ordered operator products, which does not allow simple numerical schemes. There are so-called Krylov subspace method like the short iterative Lanczos method, which offer a way out here, but these are pretty complex on their own.

The Chebychev method has some additional limitations. It works only if the Hamiltonian has either a real or a purely imaginary spectrum. This excludes important special cases like absorbing boundary conditions or open quantum systems. Also, the method requires the knowledge of exact lower and upper bounds to the spectrum, otherwise the solution may diverge.

Conclusion

As a general rule, you should always use a polynomial solver if your problem fits. In particular, when you need imaginary-time propagation to relax to the ground state or some thermal state, Chebychev polynomials are the way to go. For the many cases where polynomial solvers cannot be used, you can try out the split operator method, which tends to be surprisingly efficient, or one of the ODE solvers, which are more convenient to use. For small systems in particular, the convenience of an adaptive Runge-Kutta method easily beats the rather poor performance.

In any case, you should generally truncate your Hamiltonian. Truncation makes the propagation more stable and allows larger time steps. Plus, in many cases you only distort energy eigenstates that are poorly represented on the grid anyway.

Convergence #1: Equally-spaced grids

2020-11-24T21:35:39.361000Z

I recently got an email from a user that had tried out Wavepacket with some highly accurate data, and found that the results did not match the data. While answering the email, my answer turned out to be a bit long. And I figured out that I may have accumulated quite some experience with respect to converging (wavepacket) calculations over the years. Also, these sort of issues come up every now and again, and there have been similar questions in the past. So I wanted to shape this knowledge into some blog posts about convergence, and how to minimise errors in a wavepacket calculation.

While writing this first blog post, it turned out that properly presenting the maths is rather involved. For this reason, in this post, I only discuss a very special problem: What happens if you have an equally-spaced grid and choose the grid parameters incorrectly? The advantage of equally-spaced grids is that there is a simple and intuitive answer to this question. In later posts I would try to consider other grids, where things become a lot less clear, and also consider time-step problems.

Equally-spaced grids

An equally-spaced grid is determined by three parameters: The center of the grid x_0, the interval length L, and the number of points N.
The latter parameter can also be expressed by the grid spacing . The grid points are equally distributed in the interval , including the two end points.

If you sample any wave function on these grid points, you can do a discrete Fourier transformation to get Fourier components also on an equally-spaced grid in the momentum space. The wave vectors range from -k_max to +k_max with (details depend on whether N is even or odd). For this reason, any wave function that is represented exactly at the grid points (the discrete variable representation DVR) can be equivalently written as a sum of discrete Fourier components (the finite basis representation FBR). For more details, see the wiki page on DVR methods.

In particular, it is common nowadays to use the Fourier method in wavepacket calculations. Here, the momentum operator is applied by transforming into the Fourier space with a discrete Fourier transformation, multiplying each Fourier component with the appropriate momentum and transforming back. Using the Fourier method automatically implies periodic boundary conditions, since the plane waves are generally not limited to the interval L.

Now if we assume for a moment that you manage to choose the parameter x_0 wisely, we are left with the two remaining parameters L and N, which lead to different convergence problems.

Effect of too small interval

If the interval L is chosen too narrow for the dynamics that take place, your wavepacket may travel outside of your grid. When this happens, due to the periodic boundary conditions, the wavepacket may enter the grid from the other side.

For a simple example, let us take a free particle, give it some finite momentum and propagate. The Python code using Wavepacket can be found here, and the resulting video can be seen here.

The initial state is a Gaussian with some positive momentum. As expected, it initially broadens and moves with positive velocity. However, as soon as it reaches the grid boundaries, we can see it entering the grid again thanks to the periodic boundary conditions. Also, towards the end of the video, the Gaussian has broadened so much that the right side (roughly corresponding to the components with the larger momentum) overtakes and interferes with the left side (roughly corresponding to components with smaller momentum).
The resulting oscillations are also typical artefacts if the grid is chosen too small.

In practice, you often encounter another pattern. As an example, let us dissociate a molecule with a Morse potential parametrized for an OH bond. To have a clear wavepacket, the molecule gets a pretty strong kick initially. Again, you can find the code and the resulting video.

What we see here is an initial state with enough energy to dissociate (i.e., a wavepacket travelling to infinity). However, as it reaches the end of the grid, due to the periodic boundary conditions, the wavepacket suddenly feels the influence of the Coulomb barrier at small internuclear distances. As a result, the wavepacket does not reenter the grid from the other side, but is reflected at the boundary of the grid. The barrier here is rather steep, which is a feature not well supported by an equally-spaced grid, but that we can discuss further below.

A problem that you may sometimes encounter is that your grid is always too small. Take the Morse oscillator example, and imagine a laser-driven dissociation of a molecule. The dissociation may take some time, so you may need to propagate for large times. But in these large times, the parts of your wavepacket that dissociated first may reach the boundary of your grid, no matter how large you choose it. What to do then?

Absorbing boundary conditions to the rescue

The solution for this problem are absorbing boundary conditions, typically taking the form of negative imaginary potentials (NIPs). These absorb the wavepacket at the locations where they are nonzero, before the wavepacket reaches the grid boundaries. There are, however, two major drawbacks:

If the potential is too large, the NIP itself can modify the dynamics.
Consider a case where the NIP at some grid point absorbs effectively the whole wavepacket in a small timestep. This effect is similar to an infinite potential wall, which would also lead to a density of zero at the grid point, and may thus reflect the wavepacket. A consequence is that you should choose the NIP to be as small and smooth as possible, and you may need to converge this as well.
With the NIP term, the Hamiltonian is no longer self-adjoint.
This has serious consequences, because some propagation schemes may no longer work. In particular, the expansion in Chebychev polynomials is only correct for a self-adjoint operator, otherwise you loose all the guarantees (error bounds) that this expansion offers. The workaround here is to move the NIP out of the Hamiltonian, and absorb the wavepacket after every propagation step. This is currently not implemented in the C++ version of Wavepacket, though, because it requires exponentiation of operators.

Let us take the Morse oscillator from the previous section as an example. We might say that an OH distance of 8 atomic units (4 Angstroem) is almost dissociated. Also, from studying the dynamics, it takes at least 100 atomic units (2.5 fs) for the shoulder to travel from 8 a.u. to about the end of the grid. For a constant potential, the absorption of the density is roughly (the factor two comes because the absorption acts on the wavefunction, but the density is the square of the wavefunction). If we want to absorb 99% of the density (A = 0.01), we get a value of V_0 approximately 0.025 H. In practice, we might want to use a smooth turn-on (e.g., harmonic potential) that is zero at r = 8 a.u, and about the intended value at r = 10 a.u., that is, in atomic units, .

Applying this absorbing potential to the Morse oscillator demo gives the results shown in this video (the Python script). For clarity, the y-range has been reduced considerably. We see how the wavepacket is heavily absorbed up to about two percent of the total density. You might want to increase the grid or play around with the NIP if this is not enough, but the principle works as the back of the envelope calculation suggested.

Effect of too few grid points

For a given interval length, the equally-spaced grid also describes an equally-spaced grid of points in the Fourier space. The grid in the Fourier space spans roughly the interval [-k_max, k_max] with . (Details differ for example between even and odd N.)

If this grid is not sufficient to contain the wavepacket, you get artefacts known as aliasing. Roughly, the plane wave and give exactly the same function values at the grid in the real space. But the discrete Fourier transformation always maps this oscillation to the lower wave number k. As a consequence, we have the same effect as the real-space periodic boundary conditions: If the wavepacket accumulates too large Fourier components k > k_max (k < -k_max), it "reenters" the grid at the other side of the Fourier plane

As a simple example, let us consider a linear ramp. A Gaussian wavepacket on the linear ramp experiences some broadening, but also accelerates according to classical equations of motion. If the ramp is sufficiently steep, the wavepacket accumulates sufficient momentum to exceed the range in the Fourier plane.

Let us look at an example, see the code here, and the realspace video here. We can see how the initial wavepacket accelerates towards the negative end of the grid. Shortly before reaching the end of the grid, however, something happens and the wavepacket travels back (upwards on the potential ramp!). The effect is easily understood if we look at the wavepacket in the Fourier space. We see how the wavepacket accumulates momentum until it reaches the end of the grid in Fourier space, and reenters on the side with the positive momenta. After the whole Gaussian has traversed this boundary, we have a Gaussian that has a positive momentum and travels the ramp upwards. This slowly decreases the momentum, and the whole process starts anew.

Hence, equally-spaced grids provide a very simple measure for convergence: eyeballing. Always make sure that your wavepacket stays sufficiently far away from the grid boundaries in Fourier space. Now what is "sufficiently far away"? In principle, the exact definition is a bit complicated and related to off-diagonal operator elements in Fourier space (see the explanation of the DVR approximation). However, for most practical cases, "sufficiently far away" is the typical extension of the Fourier transform of the potential.

Takeaway message

Checking the quality of an equally-spaced grid requires observation of the dynamics both in real space and in Fourier space. This should always be done. In particular, if you find indications of parts of your wavepacket hitting the grid boundaries, you need to increase the interval or number of grid points. These boundary effects are, however, easy to understand and spot. For multi-dimensional problems, you should probably monitor the reduced density.

If there are structural problems, and you cannot make the grid large enough, you can always try to use absorbing boundary conditions. In theory, you could also employ them in the Fourier space, although this seems rather uncommon.

So what is my personal takeaway?
For one thing, writing up such an article is surprisingly hard work. On the other hand, it came back to me that the C++ version of Wavepacket has very poor support for exactly this use-case. The next version (that would be 0.3.3) should definitely feature some usability improvements for plotting and convergence monitoring.

Building a Python interface to a CMake library

2020-10-08T20:42:23.007000Z

I have recently rewritten the way how the Python module for Wavepacket is built, and learned a lot about Python along the way. The learning curve was rather steep, and it felt a lot like learning the intricate details of Nuget and .Net assembly loading (this is not a compliment). However, while extension libraries for Python may not be ideally documented, the implementation is pretty robust and neat. So I thought, I'd write up what I tried and found out to assist others in a similar situation.

Side note: I do appreciate feedback.

To briefly summarize my initial situation: I have a C++ library, WavePacket, for which I want to provide an easy-to-install and easy-to-use Python interface using Pybind11.

Things start out simple ...

At the end of the day, Python's approach to loading modules, which may include libraries with native code, is pretty robust. Neglecting topics like shadowing where multiple modules have the same name, the algorithm works roughly like this:

Python has search paths where it looks for modules. You can extend these paths either by setting the environment variable PYTHONPATH to a colon-separated list of directories, or within Python by adding a path string to the list in sys.path.
When you import a module, Python looks in every search path directory. It searches for either a library with the module name ("lib<module>.so" with various variations under Unix, "<module>.pyd" under Windows) or a subdirectory with the module name.</module></module>
When it finds a library, it looks for an exported symbol "PyInit_<module>" with a certain signature, which then tells Python about the functions, classes etc. that it exports.. Such a library can be constructed relatively easily for example with the help of pybind11.</module>
If a subdirectory with the name of the module exists, it has to contain a file init.py, which is evaluated on loading the module.

One interesting feature of init.py is that it can load more modules, also relative to the module's path.

...but they can get pretty complicated

On top of all this functionality, however, there is the Python package management system (pip), which provides a simple-to-use repository of common
packages and so on. This part I was actually much less interested in, but it was way easier to find documentation on how to create a Python package then how the whole thing works. So my original approach was to bundle the Python interface in a package.

Python used to have a big blob called setuptools that would magically do all the building, packaging, installation, distribution and so on for you. Using internet sources, I finally cobbled together the following solution:

You build and install Wavepacket using CMake with whatever backend (e.g., make).
As a side effect of the build, my CMake scripts output a setup.py file for the Python compilation in the build directory.
You run this setup.py script to build and install the Python bindings
Afterwards, you can run the various Python integration tests in the build directory using CMake tools again (ctest).

The good thing was that this solution worked. However, it had some glaring deficits. On one hand, if you wanted to test that your Python interface worked, you had to use CMake/make, then Python setuptools, then CMake again, which is ugly. What was worse, especially when compiling under Windows, there were two different build systems (CMake, setuptools) with different compiler flags and potentially compilers and downright different logic.

So I decided to rewrite this part.

Scikit

Fortunately, since my first try, the Python world has moved on. Basically, the moloch setuptools has withdrawn to provide only framework functionality for the packaging. Other tools have sprung up to provide the actual frontend that deals with, for example, the compilation of a Python extension library.

An interesting approach is used by scikit-build. The basic idea here is that you already have a library that is built with CMake, and add the Python bindings there. Scikit-build then drives the CMake build, which includes the Python bindings, takes the artefacts that come out, and wraps them in a Python package. While this sounds already great, scikit-build goes further. It provides some CMake extensions (/functions) for building and linking against a given Python
installation, and, when calling CMake, sets the module path so that these extensions are easily found.

With such nice tooling at my side, my idea was to remove the split between the Python bindings and the C++ library. If the CMake build was driven by scikit-build, I would directly compile the Python bindings together with the actual wavepacket code into a single library, and have it packaged and installed. So I would have a Python setup script on top of my CMake build that would drive the compilation if you want to get out a Python package.

Alas, while trying, I discovered some stepping stones that eventually made me give up this approach:

First of all, the documentation does not cover all the functionality, so when you have some edge cases, you are left with trying out things or reading the code.
One of the issues I learned by trial and error is that scikit-build expectes a very particular directory layout. When you want to install a module named "wavepacket", you must have certain data (such as an init.py file) in a directory called "wavepacket", otherwise things did not turn out well by default. In other words, follow the published examples to the letter.
Another issue were the CMake extensions themselves. To use the full functionality, I had to compile the Wavepacket library as a module (CMake terminology for a library that is loaded at run-time). But CMake correctly prohibits linking against modules, so for example all my C++ unit tests would have to be disabled as well when compiling for Python. Also, CMake supports two different notations for declaring link dependencies, but you may only pick one per target. Scikit-build uses the legacy notation, but I prefer the modern one for boost dependencies, which also clashes then.
setuptools originally had some functionality to run unit tests, which is not present in scikit-build. While I fully understand and support the motivation (different concerns should be addressed by different tools), I need yet another tool (e.g., tox) besides CMake and scikit-build for driving my Python tests.

So while there was no real unsurmountable obstacle, things became annoying enough that I was looking for alternatives. Mind you, if you keep these things in mind, scikit-build does look like a viable alternative.

However, while playing around with scikit-build, I eventually understood enough of Python's internals to be able to build my own module by hand. Furthermore, pondering the problem, it occured to me that I do not want to provide a Python package at all. Wavepacket has some pretty non-standard dependencies (namely the underlying tensor library) that are not easily included in a distribution. So it seems more honest not to build a Python package, but only a module that can be loaded from Python.

The final approach

And that is where my final concept ended up: I use only a CMake build, which creates the Python bindings as a byproduct. The drawback is that this can
become tricky for edge cases (multiple Python installations and such), but CMake has facilities to help there.

As a quick run-down: My top-level CMakeLists.txt offers an option to build the Python bindings and searches for Python

# ...
option(WP_BUILD_PYTHON "Build the Python interface" ON)
# ...
if (WP_BUILD_PYTHON)
     find_package(Python3 COMPONENTS Interpreter Development)
     if (NOT Python3_FOUND)
          message(ERROR "Need Python3 interpreter and development environment to compile Python module")
     endif()
endif()
#...
add_subdirectory(src)
if (WP_BUILD_PYTHON)
     add_subdirectory(python)
endif()

It then descends into two subdirectories: src, which builds the C++ library with the Wavepacket functionality, and python, which builds the Python bindings.

The CMakeLists.txt under python/ does some general setup (https://sourceforge.net/p/wavepacket/cpp/git/ci/master/tree/python/CMakeLists.txt):

It checks for the existence of certain Python modules. This can be done by importing the module with a Python command line and checking the return value
of the interpreter.
It queries the include directories for Pybind11. These are fortunately accessible from Python.
It adds a CMake test to run the various Python tests. This requires some fiddling, because we need to tell Python where to find our Python bindings once they are compiled, and I did not quite grasp the idea of Python's unittest package, so I execute it in the directory of the tests. Also, I copy the tests to the binary directory where everything is built to to avoid cluttering my source directory with residues.

Finally, under python/ I create another directory wavepacket/ that will hold my Python bindings. The directory contains the source code for the Python bindings, a CMakeLists.txt to drive the compilation and installation, and an init.py script.

The compilation of the Python bindings is pretty standard now. The rest is a bit particular:

configure_file (
     "${CMAKE_CURRENT_SOURCE_DIR}/__init__.py"
     "${CMAKE_CURRENT_BINARY_DIR}/__init__.py"
     COPYONLY)

# Installation rules
# Note that we need to link against the wavepacket library, so we need to set an rpath
set(pythonInstallDir ${CMAKE_INSTALL_LIBDIR}/wavepacket_python/wavepacket)
set_target_properties(wavepacket_python PROPERTIES INSTALL_RPATH $ORIGIN/../..)

install(TARGETS wavepacket_python
     LIBRARY DESTINATION ${pythonInstallDir})
install(FILES __init__.py
     DESTINATION ${pythonInstallDir})

The first configure_file copies the init.py to the binary directory. As a consequence, once I build the Python bindings, I only need to add
${CMAKE_BINARY_DIR}/python to PYTHONPATH and the Python installations can find my wavepacket package just as if it had been installed. This is useful for running the Python tests and demos directly in the binary directory before installation.

When installing the libraries under ${Install}/lib, I decided to put the Python files under ${Install}/lib/wavepacket_python/wavepacket. This way, you always add ${Install}/lib/wavepacket_python to your Python path and can access the "wavepacket" module. Into this directory, I install my Python binding library
and the init.py file.

One last thing to be taken into account is that my Python binding library needs to know the location of the C++ library. This I solved by adding a relative ("origin") rpath that tells the loader where to look for other libraries; this feature should be supported on most modern Unix variants. A Windows build would need some special processing, but building Wavepacket under Windows is nothing a normal user would want anyway.

Now the init.py file is pretty straight-forward:

from .wavepacket_python import *

It just directs Python to the wavepacket_python module (i.e., library) in the current directory (hence the dot) and imports all symbols from there. This nicely decouples the name of the library with the Python bindings (libwavepacket_python.so) from the name of the Python module (encoded in the directory name).

While this solution is still quite some way from perfect, it solves my initial problems. I can now compile and isntall the C++ library and the Python bindings with a single CMake call, and also run all the tests using only CMake.

WavePacket 0.2.4 released

2019-04-28T15:33:49.943000Z

I just released the next C++-version of WavePacket. Some excerpts from the NEWS file:

renamed some key classes to disentangle the concept "operator" from "expression/equation/Liouvillian". The new naming system is hopefully much more intuitive. See this blog post for the explanations.
The Python interface can be installed as a regular Python module. There are still some loose ends here and there, but it works.
added a couple new demos
- highlighted demo that explains how to setup thermal states (propagation in imaginary time, random wave functions etc.). See here.
- demo that uses multiple electronic states (MolElectronic/OH/1 and 2)
- unfinished demo that sets up a Lindblad relaxation (MolVibration/OH/3)
improved several PropagationObservers
- LoggingObserver deals correctly with multiple electronic states
- Plot1DObserver is way easier to setup, and can plot output for multiple electronic states
added a projection operator
some other minor fixes and improvements

Recent posts to blog

Perils of compiling software under Windows

Introduction

Problems (and solutions) when compiling under Windows

1. No package management

2. The dynamic loader is comparably dumb

3. The dynamic loader may get in the way

4. Licenses

5. Miscellaneous other issues

Conclusions

Wavepacket 0.3.6 released

1. Windows binaries and simpler compilation

2. 2D contour plots

3. A solver based on Faber polynomials

4. Smaller features

Outlook

Wavepacket 0.3.5 released

1. Much better support for open-system dynamics

2. More consistent names

3. Build cleanups

4. Smaller features

Outook

Wavepacket 0.3.4 released

Dependency management with DotNet or Why you should use Paket

Shortest possible summary

1. Dependency management is hard

Example 1

Example 2

What to do with a diamond dependency

2. Dependency management in DotNet is harder

3. Ditch Nuget (the tool), choose Paket instead.

Wavepacket 0.3.3 released

Convergence #2: time steps

Some theory

Model system

Finding a good time step and monitoring convergence

Application to model system

Corollary

Solution schemes

Runge-Kutta and other ODE solvers

Advantages

Disadvantages

Split-operator method

Advantages

Disadvantages

Polynomial methods

Advantages

Disadvantages

Conclusion

Convergence #1: Equally-spaced grids

Equally-spaced grids

Effect of too small interval

Absorbing boundary conditions to the rescue

Effect of too few grid points

Takeaway message

Building a Python interface to a CMake library

Things start out simple ...

...but they can get pretty complicated

Scikit

The final approach

WavePacket 0.2.4 released