6
$\begingroup$

I am trying to solve for the allowed wavefunctions and energies for a 1D quartic potential well.

To do this I am using the patching method (https://engineering.dartmouth.edu/microeng/otherweb/henning/papers/jap-sub95.pdf). I integrate from x = 0 to x = patch and from x = patch to x = end with the Numerov method. I then find the zeros of the function

$$\psi_{left}(patch) *\frac{\mathrm{d} \psi_{right}(patch) }{\mathrm{d} x} - \psi_{right}(patch) *\frac{\mathrm{d} \psi_{left}(patch) }{\mathrm{d} x} $$

where left means integration from x = 0 to x = patch, and right means integration from x = end to x = patch.

I am able to find the correct wavefunctions for a quartic well of small depth. The potential:

enter image description here

And the first four wavefunctions that were found:

enter image description here

The issue I am having is that when I increase the well depth, I am unable to find the correct wavefunctions. For the following deeper well:

enter image description here

The first found wavefunction (below) for this deeper well is clearly not the first allowed energy because it has too many zero crossings:

First found wavefunction for deeper well

When I inspect the minimizing function that I find the roots of, I see that it does not cross zero at low energies.

I suspect the issue for this is that in the deeper well, the two potential wells can be considered separately because there is so little wavefunction overlap. I would, however, like my solver to be robust enough to find the isolated wavefunctions, or at the very least recognize when the two potential wells can be considered separately. I would also like my solver to solve asymmetric wells, which I am not sure it can do right now. Is there a way I can adapt my code to solve these two issues?

Please find my code attached: https://drive.google.com/drive/folders/0B1OuUJ2C7IHAcHRkcm1JLVpLMWM

Also, here is the MATLAB code reproduced below. The main function is the first code block and is the script.

%% Outputs
% E_allowed; % Allowed eigenvalues
% wavefunctions_allowed; % Allowed eigenfunctions

%% Inputs
potential_function; % potential function, energy in eV
z_axis; % the 1D problem
h; % z_step, z(2) - z(1), distance in meters
mass_particle;
E_guess_max; E_guess_min;
E_guess_step;
E_disp;

%% Set up
% Find eigenfunctions and eigenenergies for an arbitrary potential function
if(~any(strcmp(who,'MAX_SHRO_LOOP')))
    MAX_SHRO_LOOP = 200; % Maximum runs for the bisection loop
end

if(~any(strcmp(who,'TWO_PASS_PATCHING')))
    TWO_PASS_PATCHING = 1; % Maximum runs for the bisection loop
end

if(~any(strcmp(who,'E_guess_step')))
    E_guess_step = 1e-4; % MUST BE SET ACCORDINGLY
    %MUST BE SMALLER THAN THE SPACING BETWEEN ADJACENET ENERGY LEVELS
    %MUST BE LARGE SUCH THAT WE INCREMENT E AFTER EIGEN FOUND WE SEE CHANGE
end

if(~any(strcmp(who,'z_patch_vec')))
    % the x patches that can be tried
    z_patch_vec = 3 : numel(z_axis) - 2;
end

if(~any(strcmp(who,'z_patch_index')))
    z_patch_index = round(numel(z_patch_vec) / 2); % two pass z patching
end

if(~any(strcmp(who,'init_small')))
    init_small = 1e-100; % have to adjust to get good results
    % if determinant is too small then increase this number--if the eigenvalues
    % are found too quickly
    % if NaNs in wavefunctions then decrease
end


if(~any(strcmp(who,'num_to_find')))
    % adaptive eigenenergy finding...
    num_to_find = [30 60 100 500 Inf];
end

if(~any(strcmp(who,'accuracy_to_find')))
    accuracy_to_find = [1e-16 1e-8 1e-8 1e-8 1e-8]; 
    % after each num_to_find has been reached, change the accuracy
end

% E in eV
if(~any(strcmp(who,'E_guess_min')))
    E_guess_min = min(potential_function);   
end


if(~any(strcmp(who,'E_guess_max')))
    E_guess_max = max([potential_function(1) potential_function(end)]);    
end

if(~any(strcmp(who,'E_disp')))
    E_disp = 0;
end




%%
mass_particle = mass_particle / q;

Es_guessed = [];
crossfires_guessed = [];

% the calculated stuff...
E_allowed = [];
wavefunctions_allowed = [];
left_allowed_wavefunctions = [];
right_allowed_wavefunctions = [];

E_guess = E_guess_min;
E_bracket_low = E_guess_min;
E_bracket_high = E_guess_min;

% for displaying stuff
E_guess_range = abs(E_guess_max - E_guess_min);
curr_percent = 10;
fprintf('\nShrodinger: ');

z_patches = []; % the x_patches tried
if_loop_counter = 1; % the counter for the adaptive x pacthing

while E_guess < E_guess_max

    % adapative eigenenergy error finding
    EIGEN_ERR = accuracy_to_find(sum(numel(E_allowed) > num_to_find) + 1);  
    % adaptive x patching
    z_patch = z_patch_vec(z_patch_index);
    % print progress bar
    if( 100 * abs(min(potential_function) - E_guess) / E_guess_range >= curr_percent )
        fprintf(1,'%2d%%... ',(curr_percent));
        curr_percent = curr_percent + 10;
    end

    % Find a good x_patch value

    % Low energy bracket
    T = h^2 * 2 * mass_particle * (potential_function - E_bracket_low) / (12 * h_bar^2); % T function
    guess_left_bracket_low  = Numerov_Left(0, init_small, T);
    guess_right_bracket_low = Numerov_Right(0, init_small, T);
    [crossfire_mat_bracket_low val_low] = make_crossfire_mat(guess_left_bracket_low, guess_right_bracket_low, z_patch, h);

    % high energy bracket
    E_bracket_high = E_bracket_low + E_guess_step;
    T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function
    guess_left_bracket_high  = Numerov_Left(0, init_small, T);
    guess_right_bracket_high = Numerov_Right(0, init_small, T);
    crossfire_mat_bracket_high = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
    %--this code can be deleted...

    % poor patching crossfire returns NaN, try another patching value
    if(isnan(crossfire_mat_bracket_low) ||  isnan(crossfire_mat_bracket_high))

        if_loop_counter = 1 + if_loop_counter;
        z_patch_index = mod(z_patch_index + 1, numel(z_patch_vec)) + 1;
        z_patch = z_patch_vec(z_patch_index);

        if(if_loop_counter > numel(z_patch_vec))
            disp('patching'); % all the x patchings were tried, none worked...
            break;
        end

    end
    %--this code can be deleted...

    %--this if statement (but not the code in the block) can be deleted
    % only continue if the x patching is okay and neither crossfire is not NaN
    if(~isnan(crossfire_mat_bracket_low) && ~isnan(crossfire_mat_bracket_high))

        if_loop_counter = 1; % reset the x patch counter so we can do it again later if needed
        z_patches = [z_patches; z_patch]; % the x_patches we have tried

        crossfires_guessed = [crossfires_guessed; val_low];
        Es_guessed = [Es_guessed; E_bracket_low];  

        %-------------------------- bracket guesses --------------------------%
        % Want to bound the eigenenerngy between the low and high brackets
        % Now having found a good x_patch, let us continue to bracket the eigenvalue guess

        % low energy bracket is either the min energy (if 1st try) or the
        % previous eigenenergy + guess_step

        % high energy bracket search
        % determinant has to flip signs
        while(crossfire_mat_bracket_low * crossfire_mat_bracket_high > 0)

            E_bracket_high = E_bracket_high + E_guess_step;

            T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function

            guess_left_bracket_high  = Numerov_Left(0, init_small, T);
            guess_right_bracket_high = Numerov_Right(0, init_small, T);

            [crossfire_mat_bracket_high val] = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_bracket_high];  

        end

        counter = 1; % counter for the binary search
        found_flag = 0; % 1 = found an eigenenergy

        % set the first guess to be the low bracket
        E_guess_previous = E_bracket_low;
        crossfire_mat_previous = crossfire_mat_bracket_low;

        % guess
        E_guess = (E_bracket_high + E_bracket_low) / 2;

        % binary search
        % search until MATLAB resolution reached or iteration counter satisfied
        while(abs(E_guess_previous - E_guess) > EIGEN_ERR && counter < MAX_SHRO_LOOP && ~found_flag)

            if(E_guess > E_guess_max)
                found_flag = 1;
                break
            end

            % for the bisection guess
            T = h^2 * 2 * mass_particle * (potential_function - E_guess) / (12 * h_bar^2); % T function
            guess_left  = Numerov_Left(0, init_small, T);
            guess_right = Numerov_Right(0, init_small, T);        
            [crossfire_mat val] = make_crossfire_mat(guess_left, guess_right, z_patch, h);

            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_guess];  

            % E guess correct
            if(crossfire_mat == 0)

                E_allowed = [E_allowed; E_guess];
                right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
                left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];
                found_flag = 1; % will not add to the allowed stuff at the end of the bisection loop

            % guessed too low, no sign change
            elseif(crossfire_mat * crossfire_mat_previous > 0)

                E_bracket_low = E_guess;

            else % guessed too high and the sign changed

                E_bracket_high = E_guess;

            end

            % diagnostic display
            if(E_disp)
                disp(E_guess);
            end

            % go onto the next bisection loop iteration
            E_guess_previous = E_guess;
            crossfire_mat_previous = crossfire_mat; 
            E_guess = 0.5 * (E_bracket_low + E_bracket_high); % New guess

            counter = counter + 1;            

        end

        % if the binary search stopped
        % stopped because the difference between two eigenenergies was so
        % small or because the max iterations was reached...not because
        % determinant was zero
        if(~found_flag)

            E_allowed = [E_allowed; E_guess];
            right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];    
            left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];

        end

        % increment the brackets for the next run
        E_bracket_low = E_allowed(end) + E_guess_step;
        E_bracket_high = E_bracket_low + E_guess_step;

    end

end


%%
% Clean up wavefunctions
% wavefunction correction
wavefunctions_allowed = zeros(numel(z_axis), numel(E_allowed));
patching = [];
for i = 1 : numel(E_allowed)

     [wavefunctions_allowed(:,i) patch_out] = wavefunction_patch(left_allowed_wavefunctions(:,i), right_allowed_wavefunctions(:,i), z_axis, z_patches(i), E_allowed(i), potential_function);
     patching = [patching; patch_out];

end

% flag_zerocross_eigenlevel = 1; % so far ok
% % num nodes ~= the energy level
% for i = 1 : numel(E_allowed)
%    
%     % zero crossings not commensurate with eigenvalue level
%     if(numel(find_roots(wavefunctions_allowed(:,i))) ~= i - 1)
%         flag_zerocross_eigenlevel = 0;
%         break
%     end
%     
% end
% 
% if(~flag_zerocross_eigenlevel)
%     
%     %disp('ERROR');
%     % the wells are too seperate, consider them as isolated!!
%     
% end


%% Two pass patching, now we do it again...having found the optimal patching where the agreement is closest. Necessary?
% the calculated stuff...
if(TWO_PASS_PATCHING)

    E_allowed = [];
    wavefunctions_allowed = [];
    left_allowed_wavefunctions = [];
    right_allowed_wavefunctions = [];

    % E in eV
    E_guess_min = min(potential_function);
    E_guess_max = max([potential_function(1) potential_function(end)]);

    E_guess = E_guess_min;
    E_bracket_low = E_guess_min;
    E_bracket_high = E_guess_min;

    % for displaying stuff
    E_guess_range = abs(E_guess_max - E_guess_min);
    curr_percent = 10;
    fprintf('\nShrodinger: ');

    while E_guess < E_guess_max && numel(E_allowed) < numel(patching)

        % adapative eigenenergy error finding
        EIGEN_ERR = accuracy_to_find(sum(numel(E_allowed) > num_to_find) + 1);  

        % second pass patching
        z_patch = patching(numel(E_allowed) + 1);

        % print progress bar
        if( 100 * abs(min(potential_function) - E_guess) / E_guess_range >= curr_percent )
            fprintf(1,'%2d%%... ',(curr_percent));
            curr_percent = curr_percent + 10;
        end

        % Low energy bracket
        T = h^2 * 2 * mass_particle * (potential_function - E_bracket_low) / (12 * h_bar^2); % T function
        guess_left_bracket_low  = Numerov_Left(0, init_small, T);
        guess_right_bracket_low = Numerov_Right(0, init_small, T);
        [crossfire_mat_bracket_low val_low] = make_crossfire_mat(guess_left_bracket_low, guess_right_bracket_low, z_patch, h);

        % high energy bracket
        E_bracket_high = E_bracket_low + E_guess_step;
        T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function
        guess_left_bracket_high  = Numerov_Left(0, init_small, T);
        guess_right_bracket_high = Numerov_Right(0, init_small, T);
        crossfire_mat_bracket_high = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);

        crossfires_guessed = [crossfires_guessed; val_low];
        Es_guessed = [Es_guessed; E_bracket_low];  

        %-------------------------- bracket guesses --------------------------%
        % Want to bound the eigenenerngy between the low and high brackets
        % Now having found a good x_patch, let us continue to bracket the eigenvalue guess

        % low energy bracket is either the min energy (if 1st try) or the
        % previous eigenenergy + guess_step

        % high energy bracket search
        % determinant has to flip signs
        while(crossfire_mat_bracket_low * crossfire_mat_bracket_high > 0)

            E_bracket_high = E_bracket_high + E_guess_step;

            T = h^2 * 2 * mass_particle * (potential_function - E_bracket_high) / (12 * h_bar^2); % T function

            guess_left_bracket_high  = Numerov_Left(0, init_small, T);
            guess_right_bracket_high = Numerov_Right(0, init_small, T);

            [crossfire_mat_bracket_high val] = make_crossfire_mat(guess_left_bracket_high, guess_right_bracket_high, z_patch, h);
            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_bracket_high];  

        end

        counter = 1; % counter for the binary search
        found_flag = 0; % 1 = found an eigenenergy

        % set the first guess to be the low bracket
        E_guess_previous = E_bracket_low;
        crossfire_mat_previous = crossfire_mat_bracket_low;

        % guess
        E_guess = (E_bracket_high + E_bracket_low) / 2;

        % binary search
        % search until MATLAB resolution reached or iteration counter satisfied
        while(abs(E_guess_previous - E_guess) > EIGEN_ERR && counter < MAX_SHRO_LOOP && ~found_flag)

            if(E_guess > E_guess_max)
                found_flag = 1;
                break
            end

            % for the bisection guess
            T = h^2 * 2 * mass_particle * (potential_function - E_guess) / (12 * h_bar^2); % T function
            guess_left  = Numerov_Left(0, init_small, T);
            guess_right = Numerov_Right(0, init_small, T);        
            [crossfire_mat val] = make_crossfire_mat(guess_left, guess_right, z_patch, h);

            crossfires_guessed = [crossfires_guessed; val];
            Es_guessed = [Es_guessed; E_guess];  

                % E guess correct
            if(crossfire_mat == 0)

                E_allowed = [E_allowed; E_guess];
                right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];
                left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];
                found_flag = 1; % will not add to the allowed stuff at the end of the bisection loop

                % guessed too low, no sign change
            elseif(crossfire_mat * crossfire_mat_previous > 0)

                E_bracket_low = E_guess;

            else % guessed too high and the sign changed

                E_bracket_high = E_guess;

            end

            % diagnostic display
            if(E_disp)
                disp(E_guess);
            end

            % go onto the next bisection loop iteration
            E_guess_previous = E_guess;
            crossfire_mat_previous = crossfire_mat; 
            E_guess = 0.5 * (E_bracket_low + E_bracket_high); % New guess

            counter = counter + 1;            

        end

        % if the binary search stopped
        % stopped because the difference between two eigenenergies was so
        % small or because the max iterations was reached...not because
        % determinant was zero
        if(~found_flag)

            E_allowed = [E_allowed; E_guess];
            right_allowed_wavefunctions = [right_allowed_wavefunctions, guess_right];    
            left_allowed_wavefunctions = [left_allowed_wavefunctions, guess_left];

        end

        % increment the brackets for the next run
        E_bracket_low = E_allowed(end) + E_guess_step;
        E_bracket_high = E_bracket_low + E_guess_step;

    end

    %%
    % Clean up wavefunctions
    % wavefunction correction
    wavefunctions_allowed = zeros(numel(z_axis), numel(E_allowed));
    patching_second_pass = [];
    for i = 1 : numel(E_allowed)

         [wavefunctions_allowed(:,i) patch_out] = wavefunction_patch(left_allowed_wavefunctions(:,i), right_allowed_wavefunctions(:,i), z_axis, patching(i), E_allowed(i), potential_function);
         patching_second_pass = [patching_second_pass; patch_out];

    end

end

clearvars MAX_SHRO_LOOP TWO_PASS_PATCHING E_guess_step z_patch_vec z_patch_index init_small...
    accuracy_to_find num_to_find E_guess_min E_guess_max

function u = Numerov_Left(init1, init2, T)
    % init1 = x(1) wavefunction value, must be zero for bound particle
    % init2 = x(2) wavefunction value, can be anything b/c normalization
    % T function for Numerov

    u = zeros(numel(T), 1);
        u(1) = init1; % boundary condition
        u(2) = init2; % boundary condition

    for i = 2 : numel(T) - 1
        % Solve for u(x + h)
        u(i + 1) = (( (2 + 10 * T(i)) * u(i) - (1 - T(i - 1)) * u(i - 1) ) / (1 - T(i + 1)));

    end

end

function [sign_val, val] = make_crossfire_mat(left, right, patch_index, h)

    crossfire = [left(patch_index) right(patch_index);...
        (left(patch_index) - left(patch_index - 1)) / h ...
        (right(patch_index) - right(patch_index -1)) / h];
    val = det(crossfire);

    sign_val = sign(val);

end

function u = Numerov_Right(init1, init2, T)
    % init1 = x(1) wavefunction value, must be zero for bound particle
    % init2 = x(2) wavefunction value, can be anything b/c normalization
    % T function for Numerov

    u = zeros(numel(T), 1);
        u(end) = init1; % boundary condition
        u(end - 1) = init2; % boundary condition

    for i = numel(T) - 1 : -1 : 2
        % Solve for u(x - h)
        u(i - 1) = (( (2 + 10 * T(i)) * u(i) - (1 - T(i + 1)) * u(i + 1) ) / (1 - T(i - 1)));

    end

end

function [out patch] = wavefunction_patch(left, right, z_axis, z_patch, eigenenergy, band_function)
% patch the left and right wavefunctions together

    % for symmemtric potentials
    if(sum(band_function - fliplr(flipud(band_function))) == 0)
        patch = round(numel(band_function) / 2);
    else

        % patch at location that is closest and nonzero
        % the division will make it nonzero because zero will go to inf!!
        left = left * right(z_patch) / left(z_patch);
        [val patch] = min((abs(left.^2 - right.^2)) ./ (left.^2));

        % patch at a peak
        temp = find_roots(gradient(left));
        patch = temp(1); % return the patching point so we can go back and rerun the crossfire with this new patching point
    end 

    out = [left(1 : patch); right(patch + 1 : end)];
    out = out / sqrt(trapz(z_axis, out.^2));

end
$\endgroup$
  • $\begingroup$ Do you need to use Numerov's method? $\endgroup$ – nicoguaro Aug 5 '16 at 2:36
  • $\begingroup$ No, do you think an alternate integration scheme would be better? $\endgroup$ – Geoffrey Xiao Aug 5 '16 at 4:56
  • $\begingroup$ It would be a good idea if you write the potential for your equation and the figures of your eigenvalues. I don't know about this method, that is why I asked. But if you can use other methods like Finite Differences, Finite Elements or Ritz method. I discussed a little bit in this answer. $\endgroup$ – nicoguaro Aug 8 '16 at 20:42
  • 2
    $\begingroup$ These kinds of questions are usually much easier to answer if you post your code directly (if possible), not just a description of the method in a paper (which isn't usually sufficient). Is it possible to see the code? $\endgroup$ – Kirill Aug 9 '16 at 1:20
  • 3
    $\begingroup$ A standard method is to use basis set expansion. Problems of the type you describe can be solved easily by using particle in a box wavefunctions as a basis, see J. Chem. Education. 2011, vol 88, p 929 for examples and sone code. $\endgroup$ – porphyrin Aug 9 '16 at 19:14
2
$\begingroup$

While I cannot help you with your specific implementation, I want to point out to an alternative method (as already indicated in a comment to phil's answer) : Marston's "Fourier Grid Hamiltonian" (FGH).

The FGH is a pseudo-spectral method and yields a simple recipe for constructing the discretized Hamiltonian matrix for bound systems. As usual, eigenvectors and eigenvalues of the discrete Hamiltonian represent discrete wavefunctions and corresponding energies. For the FGH, the Hamiltonian $H = T+V$ is split in a kinetic energy part $T$ and the potential part $V$. In position space the matrix elements $\left< x\vert V\vert x'\right>$ are just $V(x)\delta(x- x')$, so the potential energy contributes only the Hamiltonian matrix' diagonal. For the matrix elements of the kinetic energy operator, a band-limited (discrete $x_i$!) plane wave basis is used and the matrix elements $\left< x_i\vert T\vert x_j\right>$ are computed analytically.

You'll find a great discussion of this and various related Fourier methods in David Tannor's book "Introduction to Quantum Mechanics: a time-dependent perspective". Original references are also given there.

You'll find a Jupyter notebook with a simple Python/Cython implementation that tackles a quartic potential here. To change the potential you work with, find the line with pot = lambda x: x**4 - 20*x**2 and change it to any other potential you're interested in. Make sure to compute only bound states for which $\left<\psi\vert\psi\right>$ has decayed to zero close to the boundary of your $x$ domain.

$\endgroup$
0
$\begingroup$

The standard way to find the eigenvalues of the Schrodinger equation is called "imaginary time propagation". You change the coordinates, t=-i\tau, and integrate in the \tau direction. Any random initial condition will converge to the lowest energy eigenstate. The resulting equation is solved by splitting methods: First propagate the kinetic energy using Fast Fourier Transforms and then propagate the potential. This is possible since the solution decays fast (exponentially) as |x| grows and you can effectively replace your boundary conditions by periodic ones.

Here is some matlab code for a second order method

% in order to set up the problem, you need a spatial grid and a momentum grid, 
% since you will be using FFTs, the number of grid points should be a power of 2
Dx  = 0.01; 
x   = [-1:Dx:1-Dx];

% Now you define the eigenvalues in Fourier space
k           = zeros(1,Len).';
k(2:Len/2)  = 1:Len/2-1;                %  (x, 1, ..., Len/2-1, x, x, ..., x)
k((Len/2+1):Len) = -k((Len/2+1):-1:2);  %  (x, ...,x, -Len/2+1, ..., -1)
k(Len/2+1) = 0;

% from this you derive the momentum operator in Fourier space
p = -1i k.^2;

% This allows you to define the full problem, kinetic energy and potential, adjust as you please
E_kin = p.^2 / 2/mass;
V     = x.^4 - x.^2 ;

% you want to solve e^(h(T+V))Psi
% you approximate by the Strang splitting e^(h/2 V)e^(h T)e^(h/2 V)
for j=1:1:n
        Psi = ifft( Psi);
        Psi =  exp(- 1i*h* E_Kin ) .* Psi;
        Psi =  fft( Psi);
    % Apply potential V
        Psi = exp(  -1i*h/2* V)
    % Normalize the result in each step
         Psi = Psi / ( norm(Psi,2)^2*Dx);
   end
end

More sophisticated (higher order) methods can be found here, as well as many references to the problem: http://personales.upv.es/serblaza/2013JChemPhys.pdf

Since you want more eigenstates, there are several procedures: 1: iteration. compute the first eigenstate. Then, you restart the procedure, but in every step, you subtract the projection to this state. E.g. for the second state psi_2, in each step you write psi_2 = psi_2 - *psi_1.

As a complete alternative, you can discretize the Laplacian using finite differences (or use again FFTs) and then use the inverse power method. given your Hamiltonian H, iterate over n: (H-\lambda^(-1)I)^(-1)\psi_(n+1) = \psi_n for some guess of the eigenvalue \lambda. You can find details on wikipedia.

$\endgroup$
  • $\begingroup$ Are you one of the author in the reference in the hyperlink? If so, you should explicitly state it. $\endgroup$ – nicoguaro Aug 9 '16 at 18:29
  • $\begingroup$ why would that be necessary? $\endgroup$ – phil Aug 10 '16 at 2:02
  • $\begingroup$ Is to avoid self promotion. That way people know that you are the author. $\endgroup$ – nicoguaro Aug 10 '16 at 2:21
  • $\begingroup$ While imaginary time propagation is certainly an appealing method, I wouldn't call it the "standard way" of solving Schrödinger's equation. You may check e.g. Marston's Fourier Grid Hamiltonian Method (scitation.aip.org/content/aip/journal/jcp/91/6/10.1063/1.456888 or www.chem.yale.edu/~batista/classes/v572/fourier_grid_hamiltonian.pdf) for a simple yet effective&efficient method. $\endgroup$ – AlexE Aug 10 '16 at 18:38
  • $\begingroup$ Maybe I should have specified that it is the standard method if the first few eigenvalue/eigenvector pairs are sought. The method referenced by AlexE solves the entire fully discretised Hamiltonian and makes sense only if many of such pairs are needed which didn't seem to be the case here (or possibly ever since there are infinitely many anyway). $\endgroup$ – phil Aug 11 '16 at 0:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.