25

The Schroedinger equation is effectively a reaction-diffusion equation $$ i\frac{\partial\psi}{\partial t}=-\nabla^2\psi+V\psi\tag{1} $$ (all constants are 1). When it comes to any partial differential equation, there's two ways to solve it: Implicit method (adv: large time steps & unconditionally stable, disadv: requires matrix solver that can give bad ...


15

Finding the eigenvalues for the Schrödinger equation is really similar to finding the eigenvalues for the wave equation. You start with your differential equation $$\left[-\frac{1}{2}\nabla'^2 + V(r)\right]\psi(\mathbf{r}) = E' \psi(\mathbf{r})$$ where we did the change of variable $(x,y,z) \rightarrow (a_0 x, a_0 y, a_0 z)$, with $a_0 \equiv 1$ Bohr, $E'= E/...


13

The best approach is to use an ODE solver that is guaranteed to conserve the norm of the initial condition, i.e., for which $\|y_n\| = \|y_0\|$ for all $n\in\mathbb{N}$. Such solvers exist, and are called geometric integrators, since they preserve geometric properties of the exact solution (in this case, that energy is conserved, i.e., $\frac{d}{dt}\|y(t)\| =...


10

Kyle Kanos's answer looks to be very full, but I thought I'd add my own experience. The split-step Fourier method (SSFM) is extremely easy to get running and fiddle with; you can prototype it in a few lines of Mathematica and it is, extremely stable numerically. It involves imparting only unitary operators on your dataset, so it automatically conserves ...


10

If your discrete equations represent your problem properly you should have the right eigenvalues, although numerically you can get a small imaginary part when solving the eigenvalue problem. As an example, let us consider the axisymmetric case of a cylindrical well. The Schrödinger equation would be written as $$\frac{\partial^2 \psi}{\partial r^2} + \frac{...


9

256 equations is a relatively small number. All of the usual integrators, such as those included in Matlab, Maple or Mathematica should have no real problem with equations of this size and should be able to return answers in a fraction of the time it would take an algorithm you would implement yourself, because they use sophisticated explicit/implicit and ...


8

Kirill's answer is a good method but it is hampered that the fact that you need to tune it to a specific range in energy, so it is not appropriate if you have a broad-band wavepacket you need to absorb. An alternative approach is exterior complex scaling, which is reviewed well in Infinite-range exterior complex scaling as a perfect absorber in time-...


6

The solutions for the equation are in $$\psi \in \mathbb{C}^{3M}\times\mathbb{R}^+ \enspace .$$ If the number of electrons is small enough you can just use any traditional method. Like a domain discretization method (Finite Difference, Finite Element, Boundary Element), or a pseudospectral method. Since solving this equation is not more difficult than ...


6

I didn't follow exactly your notation so I can't say for sure what this would look like for your example, but two suggestions: if you really want your life made simple, check out qutip, the Quantum Toolbox in Python. It has classes specifically for quantum operators and state vectors, supports tensor product of operators in different spaces and keeps track ...


5

In the literature these boundary conditions go by the name of absorbing boundary conditions (or nonreflecting, open, radiation, invisible, far-field), and this is a well-known topic. One clear description I think is Absorbing Boundary Conditions for the Schrödinger Equation by Fevens and Jiang. Here is one approach (described in the above paper; I haven't ...


5

I can recommend using the finite-difference time-domain (FDTD) method. I even wrote a tutorial some time back that should answer most of your questions: J. R. Nagel, "A review and application of the finite-difference time-domain algorithm applied to the Schrödinger equation," ACES Journal, Vol. 24, No. 1, February 2009 I have some Matlab codes that run ...


5

Caveat: I did not read beyond the statement So, what I did was to define functions for the potential and the second derivative, and use Euler's method. So there may be other issues with your code, but this one is already fundamental. You are trying to solve a wave equation with Euler's method. This is numerically unstable. To see why, note that the ...


5

As always with numerical problems, it's a good idea to simplify things as much as possible to isolate what's really going wrong. You can simplify the problem by considering the evolution of a pure state $\psi$ instead of an entire density matrix, so solving $$i\dot\psi = H\psi$$ for the vector $\psi$ instead of the matrix $\rho$. Additionally, you can ...


5

I think that the main problem might be with the solver you are using. The Hamiltonian (matrix) in this case is Hermitian, it is even symmetric since it is purely real. You could use eigh instead of eig to take advantage of this. Furthermore, you are not removing only the first and last points but intervals of size 1 at each end. Following, I show you a ...


4

From the comments, it sounds like your problem is stiff, and using an implicit integrator will help you a lot more than trying to parallelize. But in case someone comes here looking for information on parallel time integration, you can find a discussion of some simple parallel extrapolation and deferred correction methods in this preprint of mine. It's not ...


4

I would disagree with your statement that Euler and Runge-Kutta have to be sequential. I know that you are saying this because you cannot parallelize across time steps (or across inner time steps in the case of something like RK4), however both can be parallel as long as you do so within a single time step. You would simply evaluate the derivatives of your ...


4

Numerically the sum isn't terribly easy: it behaves as $$ \sum_{n\geq 2, k\geq 1} \frac{1}{n(n^2-1)} \frac{e^{-4/k}}{(k^2-n^{-2})^2}, $$ so it converges linearly. Most numerical methods that can be expected to do well only really work for either alternating, or quickly converging sums. With a sum that converges linearly, you can only really get anywhere if ...


4

I know this problem well from my own research: it is given by the fact that the equation is very stiff. Thus, it is likely that you're doing nothing wrong (--although I haven't inspected your code). So where does the stiffness come from? The problem In order to solve the equation numerically, it is common to pick a grid and use finite-differences to ...


4

The SciPy tutorial explicitly states Note that ARPACK is generally better at finding extremal eigenvalues: that is, eigenvalues with large magnitudes. In particular, using which = 'SM' may lead to slow execution time and/or anomalous results. A better approach is to use shift-invert mode. and goes on to describe that. Basically, if $\lambda$ is the ...


4

I don't have enough reputation to comment, but I wanted to give a basic complement of gIS's answer. Simply put, if we have the decomposition over $V=V_1\otimes V_2$, which we represent with, say, the reshaping [2,4,2,4] why do we einsum over the first and third indices for trace over $V_1$ and second and fourth for trace over $V_2$? Maybe this is obvious ...


4

I think that your potential for the numerical case is wrong. The potential should be a big positive number, so the solution tends to zero outside the well when the value of the potential increases. If you check the solutions in Wikipedia for the first three (bounded) states they should be $$E_n = \frac{2\hbar^2 v_n^2}{m L^2}\, ,$$ with $v_1 = 1.28$, $v_2 =...


4

Regarding the boundary conditions: Don't be fooled by Wikipedia. Yes, the scenario in the picture suggests an absorption at the boundaries, and yes, one could use absorbing boundary conditions in order to reproduce that numerically. In the simple case of a wavepacket these are readily available, because in the end, there exists an analytical solution for the ...


4

Your problem was the lower limit of integration. It should have been $-x_e$ instead of 0, since $x_e$ is the equilibrium point for the potential and not the minimum distance. After correcting that, you get the following #%% Solution xe, lam = 1.0, 6.0 # parameters for potential xmax = 10 # Bval, Bval2 = wavefunction values at x = bound1, bound2 bound1, ...


4

When modelling spin in the Schrödinger equation, one has several alternatives which need to be chosen in advance. I'll copy an excerpt from a work of mine to give an overview (it considers atomic orbitals, but the idea is the same for a quantum well): For each alternative, the sought two-particle wavefunction would look like $$ \Psi(x_1,x_2) = \sum_{ij} c_{...


3

It is not that bosons are "encouraged" to occupy the same place, that is quite wrong. Here is an article on the spin-statistics theorem which describes the relationship between particle statistics and the transformation of the wave function under permutations of identical particles. This kind of thing requires more knowledge of quantum mechanics, so I ...


3

You want to solve for 3 to 10 particle systems (3D per particle)? As far as I am aware, mean field theories do not work especially well for so few particles, but it seems there has been DFT work on diatomic molecules. Is this a system where Born-Oppenheimer is valid? If so, I might be inclined to expand the electronic wavefunction using a linear combination ...


3

If you imagine the equations of just two bodies that are bound to each other along a fixed connection and float through space without external forces, and if you denote their centers as $x_1(t),x_2(t)$ both of which are 3-dimensional coordinates, then you can write the equations of motion as $$ \ddot x_1(t) = \lambda(x_2-x_1), \\ \ddot x_2(t) = -\...


Only top voted, non community-wiki answers of a minimum length are eligible