Integrating radial Schrodinger equation with Lennard-Jones potential using Runge-Kutta with adaptive step size ends up with a step-size of zero

Question

I'm currently taking a course in computational physics. I'm new to computational physics and programming in general. I'm using numerical recipes to try and integrate the radial Schrodinger equation with a Lennard-jones potential.

$$\left[ \frac{\hbar^2}{2m}\frac{d^2}{dr^2} + \left( E-V(r)-\frac{\hbar^2 l (l+1)}{2mr^2}\right) \right] u_l(r)=0$$

$$V(r)= \epsilon \left[ \left( \frac{\rho}{r} \right)^{12}-2\left(\frac{\rho}{r}\right)^6 \right]$$

Numerical recipes has a function called odeint which will use a fifth-order Runge-Kutta algorithm to integrate an ordinary differential equation for you. The function has an adjustable step-size which appears to be causing problems in my code. Namely, my step-size is going to zero, which causes numerical recipes to throw an error and exit prematurely.

I am doing the integration from $r_{min}=\frac{\rho}{2}$ to $r_{max}=5\rho$ numerically, and up to my minimum value analytically in order to take care of the singularity at zero. I have attached my code below, and more information on odeint and how it works can be found at: http://www.itp.uni-hannover.de/Lehre/cip/odeint_c.pdf. Can anyone help me understand where I'm going wrong?

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nr.h"
#include "nrutil.h"

#define N 2

float dxsav,*xp,**yp;  /* defining declarations */
int kmax,kount;

int nrhs;   /* counts function evaluations */

/* Schrodinger equation and L-J parameters */
double alpha = 6.12; 
double rho = 3.57;
double epsilon = 5.9;
double l = 1.0;
double energy = 3;

double rmin=1.785;
double rmax=17.85;


void derivs(float x,float y[],float dydx[])
{
    nrhs++;
    printf("xodeint check: x=%f\n", x);
    dydx[1] = y[2];
    dydx[2]=alpha*((l*(l+1)/(alpha*x*x))+epsilon*(pow(rho/x,12)-2*pow(rho/x,6))-energy)*y[1];
}

int main(void)
{
    int i,nbad,nok;
    float eps=1.0e-4,h1=0.1,hmin=0,x1=rmin,x2=rmax,*ystart;

    ystart=vector(1,N);
    xp=vector(1,200);
    yp=matrix(1,10,1,200);
    ystart[1]=0.93583;
    ystart[2]=0.17385;
    nrhs=0;
    kmax=100;
    dxsav=(x2-x1)/20.0;
//  printf("%f\n", h1);
    odeint(ystart,N,x1,x2,eps,h1,hmin,&nok,&nbad,derivs,rkqs);
    printf("\n%s %13s %3d\n","successful steps:"," ",nok);
    printf("%s %20s %3d\n","bad steps:"," ",nbad);
    printf("%s %9s %3d\n","function evaluations:"," ",nrhs);
    printf("\n%s %3d\n","stored intermediate values:    ",kount);
    printf("\n%8s %18s %15s\n","r","integral","x^2");
    for (i=1;i<=kount;i++)
        printf("%10.4f %16.6f %14.6f\n",xp[i],yp[1][i],xp[i]*xp[i]);
    free_matrix(yp,1,10,1,200);
    free_vector(xp,1,200);
    free_vector(ystart,1,N);
    return 0;
}
#undef NRANSI

This outputs

Numerical Recipes run-time error...
stepsize underflow in rkqs
...now exiting to system...

Answer 1

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 discretize the functions and operators. For simplicity, let's say your grid is equally spaces with a spacing of $\Delta x$. Usually, one wants to choose $\Delta x$ small in order to achieve a good approximation.

Now for an analysis of the terms in your equation: The kinetic energy term (i.e. the derivative) leads to upper energies of $\mathcal O(\frac{1}{\Delta x^2})$, as the largest frequency which is representable on the grid is $\mathcal O(\frac{1}{\Delta x})$. (I'm lazy so I write $\mathcal O$ here--you surely know the formulas). This term is contained in any Schrödinger equations and does not really pose a problem.

Now let's consider the potential functions. Those are diagonal on the grid, so a function $1/r^k$ will have a largest energy of $1/(\Delta x^k)$. In your equation, you have a maximum $k$ of $12$. It is this term which primarily kills your propagation.

Explanation

Why does this term poses a large problem? Consider your complete discretized Hamiltonian matrix $H$ (which is symmetric), and an expansion of your wavefunction in terms of your grid. The explicit solution of the Schrödinger equation is given by

$$ \psi(t) = \exp(i H t) \psi(t=0)\,. $$

Note that there is no time-ordering operator as the Hamiltonian is time-independent. Now make an eigendecomposition of your Hamiltonian, $H = U E U^\dagger$. Then, the exponential may be represented as $\exp(i H t) = U \exp(i E t) U^\dagger$. If you now have a very large eigenenergy $e_{max}$ in your diagonal matrix $E$, then you have a term like $exp(ie_{max} t)$ in your propagation.

Here enters the Shannon theorem (or alternatively the Fourier transformation), which states that in order to sample this function well, it should be $\Delta t \sim 1/e_{max}$. Otherwise, you get weird things like aliasing effects, which spoil your solution.

For the sake of convenience, one can approximate the maximal eigenenergies the complete Hamiltonian as the sum of the maximal eigenenergies of it's terms. With this, you see that the time step needs to be as small as $\Delta t \sim 1/\Delta x$. That is, if your $\Delta x$ is $10^{-2}$ (no units here, you know them), you'd need a time step of $\sim 10^{-24}$. Adaptive stepize integrators will thus lower and lower the stepsize, until they arrive their underflow constant.

This is a rather un-mathematical treatment so don't quote me on the numbers, but it intuitively explains the occuring effects.

Suggested solution

There is a simple alternative which I think is suited to your problem, which is often called the spectral method. It is applicable if your grid has a small to medium size -- say up to $N=10000$ gridpoints -- which I guess is the case here. Then, in fact, it is the method of choice for your problem.

The solution is then simply to calculate the matrix exponential $\exp(iHt)$ as sketched above, and apply it to your initial wavefunction vector. The costs are $\mathcal O(N^3)$ for the diagonalization, but you need to execute it only once. Plus, by this you get the exact numerical solution -- the solution all those ODE methods try to achieve.

The downside is that this approach is only applicable when the Hamiltonian is time-independent. Otherwise, when it changes over time, you'd need to diagonalize it often anew.

Answer 2

The first thing I'd try is an implicit integrator, to check for the possibility that your ODE is stiff. If you're using RK5 -- an explicit method -- to solve your ODE and the step size becomes very small, it could be that stability considerations are severely limiting your step size. If you want to use an ODE solver written in C, SUNDIALS or PETSc are good choices for production code. PETSc is very well supported. SUNDIALS has a mailing list that also gets a fair amount of traffic. The GNU Scientific Library also has some routines that are probably good for academic use, or prototyping.

You could also try a different RK implementation. I'm much more likely to trust code that's part of a library that people use and support than someone transcribing code from Numerical Recipes. Opinions on Numerical Recipes vary. The authors themselves say that criticism of bugs in the library is unfounded and list current bug reports (now obsolete, since they have a forum). However, most of the stories I've heard about Numerical Recipes are like this one -- not great. Your mileage may vary.

Finally, as a catchall, you could double-check (preferably test) to make sure that there are no bugs in the implementation of the right-hand side of your ODE, but this step would probably be my last step after ruling out the first two sources of problems.