!! This program is free software; you can redistribute it and/or modify !! it under the terms of the GNU General Public License as published by !! the Free Software Foundation; either version 2, or (at your option) !! any later version. !! !! This program is distributed in the hope that it will be useful, !! but WITHOUT ANY WARRANTY; without even the implied warranty of !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !! GNU General Public License for more details. !! !! You should have received a copy of the GNU General Public License !! along with this program; if not, write to the Free Software !! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! MODULE EXPO ! =========== ! ! This module contain the procedures that calculate the action of the exponential ! of the Hamiltonian (times -I*dt) on an input state. ! ! Note that there are three possibilities: ! (i) Setting mode = TAYLOR, one uses a truncated Taylor expansion to approximate ! exp(-I*dt*H) ! (ii) Setting mode = LANCZOS, one uses the Lanczos scheme implemented in the ! octopus code. ! (iii) And finally, setting mode = EXPOKIT, one uses the expokit package to ! perform the same task. This is a more sophisticated version of the previous ! option. ! ! We wanted to illustrate the two option that one has when writing software: ! to learn what algorithm/method one needs, and implement it, or to learn ! how those algorithms works, look for somebody else that has done it, and ! using it. Of course, the second option is highly recommended... !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module expo use mesh implicit none !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Only one public procedure: exponential. It is just a driver to the ! selected algorithm, which is decided with variable "mode". private public :: exponential, pot !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Set mode to the option you wnat to test. integer, parameter :: TAYLOR = 1, & LANCZOS = 2, & EXPOKIT = 3 integer, parameter :: mode = LANCZOS ! This is for internal use of the expokit package. real(8), allocatable :: pot(:, :) contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE EXPONENTIAL ! ====================== ! ! INPUT: ! v [real(8), dimension(n, n)] : the potential that defines the Hamiltonian. ! wf [complex(8), dimension(n, n)] : on input, it contains the wavefunction ! onto which the exponential of H is applied. ! dt [real(8)] : the subroutine should return exp[-I*dt*H]|wf> ! -------- ! OUTPUT: ! wf [complex(8), dimension(n, n)] : on output, it should contain the result ! ! This subroutine is merely a driver to the appropriate subroutine, depending ! on the value of the mode parameter, which should be set up above. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine exponential(v, wf, dt) implicit none real(8), intent(in) :: v(N, N) complex(8), intent(inout) :: wf(N, N) real(8), intent(in) :: dt select case(mode) case(TAYLOR); call exponential_taylor(v, wf, dt) case(LANCZOS) call exponential_lanczos(v, wf, dt) case(EXPOKIT) allocate(pot(n, n)) pot = v call exponential_expokit(wf, dt) deallocate(pot) end select end subroutine exponential !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE EXPONENTIAL_TAYLOR ! ============================= ! ! INPUT: ! v [real(8), dimension(n, n)] : the potential that defines the Hamiltonian. ! wf [complex(8), dimension(n, n)] : on input, it contains the wavefunction ! onto which the exponential of H is applied. ! dt [real(8)] : the subroutine should return exp[-I*dt*H]|wf> ! -------- ! OUTPUT: ! wf [complex(8), dimension(n, n)] : on output, it should contain the result ! ! This subroutine approximates exp(-I*dt*H)|wf> by using a truncated Taylor ! expansion of order whatever (four is a very good value, for some reasons...) !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine exponential_taylor(v, wf, dt) implicit none real(8), intent(in) :: v(N, N) complex(8), intent(inout) :: wf(N, N) real(8), intent(in) :: dt !!!!!!! MISSING CODE 8 !!!!!! MISSING CODE 8 end subroutine exponential_taylor !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE EXPONENTIAL_LANCZOS ! ============================= ! ! INPUT: ! v [real(8), dimension(n, n)] : the potential that defines the Hamiltonian. ! wf [complex(8), dimension(n, n)] : on input, it contains the wavefunction ! onto which the exponential of H is applied. ! dt [real(8)] : the subroutine should return exp[-I*dt*H]|wf> ! -------- ! OUTPUT: ! wf [complex(8), dimension(n, n)] : on output, it should contain the result ! ! This subroutine approximates exp(-I*dt*H)|wf> by using the Lanczos ! procedure, that is, projection onto a Krylov subspace. !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine exponential_lanczos(v, wf, dt) implicit none real(8), intent(in) :: v(N, N) complex(8), intent(inout) :: wf(N, N) real(8), intent(in) :: dt integer :: i integer, parameter :: korder = 20 real(8), parameter :: tol = 1.0e-5_8 complex(8), allocatable :: x(:, :, :), hm(:, :), expo(:, :), mat(:, :) real(8) :: nrm, alpha, beta, res complex(8), allocatable :: wsp(:) integer :: lwsp, ipiv(korder), iexph, ns, iflag, order lwsp = 4*korder*korder+7 allocate(wsp(lwsp)) allocate(x(n, n, korder), hm(korder, korder), expo(korder, korder), mat(korder, korder)) nrm = sqrt(zdotproduct(wf, wf)) x(:, :, 1) = wf(:, :)/nrm call zhpsi(v, x(:, :, 1), wf) alpha = zdotproduct(x(:, :, 1), wf) wf(:, :)= wf(:, :) - alpha*x(:, :, 1) hm = (0.0d0, 0.0d0) hm(1, 1) = alpha beta = sqrt(zdotproduct(wf, wf)) do i = 1, korder - 1 x(:, :, i + 1) = wf(:, :)/beta hm(i+1, i) = beta call zhpsi(v, x(:, :, i+1), wf) hm(i , i + 1) = zdotproduct(x(:,:, i) , wf) hm(i + 1, i + 1) = zdotproduct(x(:,:, i+1), wf) call zgemv('n', n**2, 2, (-1.0d0,0.0d0), x(1, 1, i), n**2, hm(i:i+1, i+1), 1, (1.0d0, 0.0d0), wf, 1) mat = (0.0d0, -1.0d0)*hm call ZGPADM(6, i+1, dt, mat, korder, wsp, lwsp, ipiv, iexph, ns, iflag) call zcopy((i+1)*(i+1), wsp(iexph), 1, expo(1:i+1, 1:i+1), 1) res = abs(beta*abs(expo(1, i+1))) beta = sqrt(zdotproduct(wf, wf)) if(beta < 1.0e-12_8) exit if(i>1 .and. res tol .and. beta > 1.0e-12_8) then write(*, '(a)') 'Warning: Lanczos exponential expansion did not converge' endif call zgemv('n', n**2, korder, cmplx(nrm, 0.0d0, 8), x, n**2, expo, 1, cmplx(0.0d0, 0.0d0, 8), wf, 1) deallocate(x, hm, expo, wsp) end subroutine exponential_lanczos !/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE EXPONENTIAL_EXPOKIT ! ============================== ! ! INPUT: ! v [real(8), dimension(n, n)] : the potential that defines the Hamiltonian. ! wf [complex(8), dimension(n, n)] : on input, it contains the wavefunction ! onto which the exponential of H is applied. ! dt [real(8)] : the subroutine should return exp[-I*dt*H]|wf> ! -------- ! OUTPUT: ! wf [complex(8), dimension(n, n)] : on output, it should contain the result ! ! This subroutine approximates exp(-I*dt*H)|wf> by using the Lanczos ! procedure, that is, projection onto a Krylov subspace. It makes use of ! the expokit package to do the real work. !*/!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine exponential_expokit(wf, dt) implicit none complex(8), intent(inout) :: wf(N, N) real(8), intent(in) :: dt integer nn, m, lwsp, liwsp double precision tol, anorm, s1, s2 complex*16, allocatable :: x(:), y(:) integer, allocatable :: iwsp(:) complex*16, allocatable :: wsp(:) integer i, j, nnz, itrace, iflag, iseed(4) complex*16 ZERO, ONE parameter( ZERO=(0.0d0,0.0d0), ONE=(1.0d0,0.0d0) ) external op nn = n**2 m = 8 ! I do not know what is optional here.... lwsp = nn*(m+2)+5*(m+2)**2+7 liwsp = m+2 allocate(x(nn), y(nn)) allocate(iwsp(liwsp), wsp(lwsp)) call zcopy(nn, wf, 1, x, 1) anorm = 1.0d0 !*--- set other input arguments. I also do not know which one is optimal here... tol = 1.0e-3_8 itrace = 0 !*--- compute w = exp(t*A)v with ZGEXPV ... call ZGEXPV( nn, m, dt, x, y, tol, anorm, & wsp, lwsp, iwsp, liwsp, op, itrace, iflag ) call zcopy(nn, y, 1, wf, 1) end subroutine exponential_expokit end module expo ! To define this routine (which just applies the Hamiltonian) is needed by expokit. ! Ideally one should modify expokit subroutine to avoid unnecessary copying... subroutine op(u, v) use mesh use expo complex(8), intent(in) :: u(*) complex(8), intent(out) :: v(*) complex(8), allocatable :: a(:, :), b(:, :) allocate(a(n, n), b(n, n)) call zcopy(n*n, u, 1, a, 1) call zhpsi(pot, a, b) call zscal(n*n, (0.0d0, -1.0d0), b, 1) call zcopy(n*n, b, 1, v, 1) deallocate(a, b) end subroutine op