A new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation

In this paper we present a new kind of discretization scheme for solving a two-dimensional time-independent Schrödinger equation. The scheme uses a symmetrical multi-point difference formula to represent the partial differentials of the two-dimensional variables, which can improve the accuracy of the numerical solutions to the order of Δx2Nq+2 when a (2Nq+1)-point formula is used for any positive integer Nq with Δx=Δy, while Nq=1 equivalent to the traditional scheme. On the other hand, the new scheme keeps the same form of the traditional matrix equation so that the standard algebraic eigenvalue algorithm with a real, symmetric, large sparse matrix is still applicable. Therefore, for the same dimension, only a little more CPU time than the traditional one should be used for diagonalizing the matrix. The numerical examples of the two-dimensional harmonic oscillator and the two-dimensional Henon-Heiles potential demonstrate that by using the new method, the error in the numerical solutions can be reduced steadily and extensively through the increase of Nq, which is more efficient than the traditional methods through the decrease of the step size.     

FFT-split-operator code for solving the Dirac equation in 2+1 dimensions

The main part of the code presented in this work represents an implementation of the split-operator method [J.A. Fleck, J.R. Morris, M.D. Feit, Appl. Phys. 10 (1976) 129-160; R. Heather, Comput. Phys. Comm. 63 (1991) 446] for calculating the time-evolution of Dirac wave functions. It allows to study the dynamics of electronic Dirac wave packets under the influence of any number of laser pulses and its interaction with any number of charged ion potentials. The initial wave function can be either a free Gaussian wave packet or an arbitrary discretized spinor function that is loaded from a file provided by the user. The latter option includes Dirac bound state wave functions. The code itself contains the necessary tools for constructing such wave functions for a single-electron ion. With the help of self-adaptive numerical grids, we are able to study the electron dynamics for various problems in 2+1 dimensions at high spatial and temporal resolutions that are otherwise unachievable.Along with the position and momentum space probability density distributions, various physical observables, such as the expectation values of position and momentum, can be recorded in a time-dependent way. The electromagnetic spectrum that is emitted by the evolving particle can also be calculated with this code. Finally, for planning and comparison purposes, both the time-evolution and the emission spectrum can also be treated in an entirely classical relativistic way.Besides the implementation of the above-mentioned algorithms, the program also contains a large C++ class library to model the geometric algebra representation of spinors that we use for representing the Dirac wave function. This is why the code is called “Dirac++”.

Program summary

Program title: Dirac++ or (abbreviated) d++Catalogue identifier: AEAS_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEAS_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 474 937No. of bytes in distributed program, including test data, etc.: 4 128 347Distribution format: tar.gzProgramming language: C++Computer: Any, but SMP systems are preferredOperating system: Linux and MacOS X are actively supported by the current version. Earlier versions were also tested successfully on IRIX and AIXNumber of processors used: Generally unlimited, but best scaling with 2-4 processors for typical problemsRAM: 160 Megabytes minimum for the examples given hereClassification: 2.7External routines: FFTW Library [3,4], Gnu Scientific Library [5], bzip2, bunzip2Nature of problem: The relativistic time evolution of wave functions according to the Dirac equation is a challenging numerical task. Especially for an electron in the presence of high intensity laser beams and/or highly charged ions, this type of problem is of considerable interest to atomic physicists.Solution method: The code employs the split-operator method [1,2], combined with fast Fourier transforms (FFT) for calculating any occurring spatial derivatives, to solve the given problem. An autocorrelation spectral method [6] is provided to generate a bound state for use as the initial wave function of further dynamical studies.Restrictions: The code in its current form is restricted to problems in two spatial dimensions. Otherwise it is only limited by CPU time and memory that one can afford to spend on a particular problem.Unusual features: The code features dynamically adapting position and momentum space grids to keep execution time and memory requirements as small as possible. It employs an object-oriented approach, and it relies on a Clifford algebra class library to represent the mathematical objects of the Dirac formalism which we employ. Besides that it includes a feature (typically called “checkpointing”) which allows the resumption of an interrupted calculation.Additional comments: Along with the program's source code, we provide several sample configuration files, a pre-calculated bound state wave function, and template files for the analysis of the results with both MatLab and Igor Pro.Running time: Running time ranges from a few minutes for simple tests up to several days, even weeks for real-world physical problems that require very large grids or very small time steps.References:

[1]

J.A. Fleck, J.R. Morris, M.D. Feit, Time-dependent propagation of high energy laser beams through the atmosphere, Appl. Phys. 10 (1976) 129-160.

[2]

R. Heather, An asymptotic wavefunction splitting procedure for propagating spatially extended wavefunctions: Application to intense field photodissociation of H⁺₂, Comput. Phys. Comm. 63 (1991) 446.

[3]

M. Frigo, S.G. Johnson, FFTW: An adaptive software architecture for the FFT, in: Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, vol. 3, IEEE, 1998, pp. 1381-1384.

[4]

M. Frigo, S.G. Johnson, The design and implementation of FFTW3, in: Proceedings of the IEEE, vol. 93, IEEE, 2005, pp. 216-231. URL: http://www.fftw.org/.

[5]

M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, M. Booth, F. Rossi, GNU Scientific Library Reference Manual, second ed., Network Theory Limited, 2006. URL: http://www.gnu.org/software/gsl/.

[6]

M.D. Feit, J.A. Fleck, A. Steiger, Solution of the Schrödinger equation by a spectral method, J. Comput. Phys. 47 (1982) 412-433.

     

Numerical solving of the vibrational time-independent Schrödinger equation in one and two dimensions using the variational method

A program package for variational solving of the time-independent Schrödinger equation (SE) in one and two dimensions is described. The first part of the the program package includes the fitting program (FIT) with which the ab initio or DFT calculated points are fitted to a computationally inexpensive functional form. Proper fitting of the potential energy surface is crucial for the quality of the results. The second part of the package consists of a program for variational solving of the SE (2DSCHRODINGER) using either a shifted Gaussian basis set or the rectangular basis set proposed by Balint-Kurti and coworkers [J. Chem. Phys. 91 (1989) 3571]. The third part of the program package consists of the calculation of the expectation values, IR and Raman spectra XPECT), and the visualization of results (PLOT). The program package is applied to study a quantum harmonic oscillator and an intramolecular, strong hydrogen bond in picolinic acid N-oxide. For the former system analytical solutions exist, while for the latter system a comparison with the experimental data is made. The advantages and disadvantages of the applied methods are discussed.     

A method for solving the molecular Schrödinger equation in Cartesian coordinates via angular momentum projection operators

A method for solving the Schrödinger equation of N-atom molecules in 3N−3 Cartesian coordinates usually defined by Jacobi vectors is presented. The separation and conservation of the total angular momentum are obtained not by transforming the Hamiltonian in internal curvilinear coordinates but instead, by keeping the Cartesian formulation of the Hamiltonian operator and projecting the initial wavefunction onto the proper irreducible representation angular momentum subspace. The increased number of degrees of freedom from 3N−6 to 3N−3, compared to previous methods for solving the Schrödinger equation, is compensated by the simplicity of the kinetic energy operator and its finite difference representations which result in sparse Hamiltonian matrices. A parallel code in Fortran 95 has been developed and tested for model potentials of harmonic oscillators. Moreover, we compare data obtained for the three-dimensional hydrogen molecule and the six-dimensional water molecule with results from the literature. The availability of large clusters of computers with hundreds of CPUs and GBytes of memory, as well as the rapid development of distributed (Grid) computing, make the proposed method, which is unequivocally highly demanding in memory and computer time, attractive for studying Quantum Molecular Dynamics.     

Accurate numerical method for the solutions of the Schrödinger equation and the radial integrals based on the CIP method

A new accurate numerical method based on the constrained interpolation profile (CIP) method to solve the Schrödinger wave equation for bound and free states in central fields and to calculate radial integrals is presented. The radial wave equation is integrated on an arbitrary grid system by the adaptive stepsize controlled Runge-Kutta method controlling the truncation errors within a prescribed accuracy. For the continuum orbitals in the highly oscillating region, the non-linear radial wave equation in the phase-amplitude representation is used. In the evaluation of the derivatives of the radial wave function, the potential energy is approximated by the CIP method. In addition, the radial integrals encountered in the computation of various atomic process are accomplished with the CIP method using the values and their analytical derivatives at the grids. This numerical procedure can be extended in a straightforward way to solve the Dirac wave equation.     

Arbitrarily precise numerical solutions of the one-dimensional Schrödinger equation

In this paper, how to overcome the barrier for a finite difference method to obtain the numerical solutions of a one-dimensional Schrödinger equation defined on the infinite integration interval accurate than the computer precision is discussed. Five numerical examples of solutions with the error less than 10⁻⁵⁰ and 10⁻³⁰ for the bound and resonant state, respectively, obtained by the Obrechkoff one-step method implemented in the multi precision mode, which include the harmonic oscillator, the Pöschl-Teller potential, the Morse potential and the Woods-Saxon potential, demonstrate that the finite difference method can yield the eigenvalues of a complex potential with an arbitrarily desired precision within a reasonable efficiency.     

Binary SCF: GAMESS improvements for energy evaluation based on SCF methods

This paper describes a new algorithm that improves computational efficiency for energy evaluation in GAMESS, named Binary SCF (BitSCF). This algorithm is based on the Conventional Self Consistent Field (CSCF) method and it presents two main innovations. The first one is the use of a bit map for the control of the 2e⁻ integrals evaluated and previously stored in disk. The bit map reduces disk demand significantly when storing 2e⁻ integrals. The second innovation is the use of i-loop balancing, a new proposal for static workload distribution through the nodes. The i-loop balancing optimizes parallelism in GAMESS, reducing an important sequential portion of its source code. The results obtained confirm a reduction of up to 49% in the disk demand for the 2e⁻ integrals storage, considering different simulated molecules on a 13-node Beowulf cluster and 6-31G and STO-6G basis set. The total and CPU times stay below the times obtained with the GAMESS CSCF method for larger molecules. The total execution and CPU times were equivalent at the CSCF times for smaller molecules. The results show an improvement in the GAMESS scalability as well as a better cost vs. benefit relationship when executing the program.     

Collision dynamics of elliptically polarized solitons in Coupled Nonlinear Schrödinger Equations

We investigate numerically the collision dynamics of elliptically polarized solitons of the System of Coupled Nonlinear Schrödinger Equations (SCNLSE) for various different initial polarizations and phases. General initial elliptic polarizations (not sech$sech$-shape) include as particular cases the circular and linear polarizations. The elliptically polarized solitons are computed by a separate numerical algorithm. We find that, depending on the initial phases of the solitons, the polarizations of the system of solitons after the collision change, even for trivial cross-modulation. This sets the limits of practical validity of the celebrated Manakov solution. For general nontrivial cross-modulation, a jump in the polarization angles of the solitons takes place after the collision (‘polarization shock’). We study in detail the effect of the initial phases of the solitons and uncover different scenarios of the quasi-particle behavior of the solution. In majority of cases the solitons survive the interaction preserving approximately their phase speeds and the main effect is the change of polarization. However, in some intervals for the initial phase difference, the interaction is ostensibly inelastic: either one of the solitons virtually disappears, or additional solitons are born after the interaction. This outlines the role of the phase, which has not been extensively investigated in the literature until now.     

Numerical errors resulting from finite-difference approximation in computations of a one-dimensional Schrödinger equation with the Trotter formula

The parallel algorithm for solving time-dependent Schrödinger equations devised by De Raedt and based on the Trotter formula is not only simple but also unconditionally stable, explicit, and local. We consider the numerical errors resulting from the finite-difference approximation of De Raedt's algorithm by comparing an exact solution of a free particle with the approximate solution calculated by using the Trotter formula, which depends on the size of the spatial-temporal lattice.     

Solving the vibrational Schrödinger equation on an arbitrary multidimensional potential energy surface by the finite element method

A computational protocol has been developed to solve the bounded vibrational Schrödinger equation for up to three coupled coordinates on any given effective potential energy surface (PES). The dynamic Wilson G-matrix is evaluated from the discrete PES calculations, allowing the PES to be parametrized in terms of any complete, minimal set of coordinates, whether orthogonal or non-orthogonal. The partial differential equation is solved using the finite element method (FEM), to take advantage of its localized basis set structure and intrinsic scalability to multiple dimensions. A mixed programming paradigm takes advantage of existing libraries for constructing the FEM basis and carrying out the linear algebra. Results are presented from a series of calculations confirming the flexibility, accuracy, and efficiency of the protocol, including tests on FHF⁻, picolinic acid N-oxide, trans-stilbene, a generalized proton transfer system, and selected model systems.     

Highly precise solutions of the one-dimensional Schrödinger equation with an arbitrary potential

This work concerns the obtaining of a highly accurate solution of the one-dimensional Schrödinger equation with a practically arbitrary potential. The approach is based on the power series method and it is implemented in the ultra high precision mode. It is shown that such an approach yields not only highly precise values of the energies but also accurate wave functions.     

Optimized implementations of rational approximations—a case study on the Voigt and complex error function

Rational functions are frequently used as efficient yet accurate numerical approximations for real and complex valued special functions. For the complex error function w(x+iy), whose real part is the Voigt function K(x,y), the rational approximation developed by Hui, Armstrong, and Wray [Rapid computation of the Voigt and complex error functions, J. Quant. Spectrosc. Radiat. Transfer 19 (1978) 509-516] is investigated. Various optimizations for the algorithm are discussed. In many applications, where these functions have to be calculated for a large x grid with constant y, an implementation using real arithmetic and factorization of invariant terms is especially efficient.     

Symplectic and multi-symplectic methods for coupled nonlinear Schrödinger equations with periodic solutions

We consider for the integration of coupled nonlinear Schrödinger equations with periodic plane wave solutions a splitting method from the class of symplectic integrators and the multi-symplectic six-point scheme which is equivalent to the Preissman scheme. The numerical experiments show that both methods preserve very well the mass, energy and momentum in long-time evolution. The local errors in the energy are computed according to the discretizations in time and space for both methods. Due to its local nature, the multi-symplectic six-point scheme preserves the local invariants more accurately than the symplectic splitting method, but the global errors for conservation laws are almost the same.     

Application of efficient composite methods for computing with certainty periodic orbits in molecular systems

Recently, we have proposed a technique for the computation of periodic orbits in molecular systems, based on the characteristic bisection method [Vrahatis et al., Comput. Phys. Commun. 138 (2001) 53]. The main advantage of the characteristic bisection method is that it converges with certainty within a given starting rectangular region. In this paper we further improve this technique by applying, on a surface of section of a Poincaré map, an iterative scheme based on the composition of the characteristic bisection method with other more rapid root-finding methods such as Newton's or Broyden's methods. Thus, the composite schemes compute rapidly with certainty periodic orbits of molecular systems. By applying these methods to the LiNC/LiCN molecular system we obtain promising results. We have reproduced previous results using considerable less CPU time. Also, we have located and computed new asymmetric families of periodic orbits.     

The MATHSCOUT Mathematica package to postprocess the output of other scientific programs

mathscout is a mathematica¹ package to postprocess the output of other programs for scientific calculations. We wrote mathscout to import data from a major program for ab initio computational chemistry into mathematica, so that we could postprocess the chemical results. It can be used to import the output of many other packages that are used, e.g. in molecular dynamics, crystallography, spectroscopic analysis, metabolic and physiological modeling, meteorology and other areas of environmental science, cosmology and particle physics. mathscout assigns a name to each table and non-tabular datum that it extracts. This name is constructed mechanically from the identifier or phrase that precedes or follows or embeds the item in the output that mathscout processes. A selection of non-contiguous items, or all the items in a section of the file, or in the entire file are extracted using simple commands. So far, we have focused on our immediate needs to postprocess the output of the Gaussian² program. Calculations on several molecules that illustrate the usage of the package are presented here and in the Supplementary Information. mathscout is shortened to msct in the software.

Program summary

Program title: msct.mCatalogue identifier: ADZQ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADZQ_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 30 396No. of bytes in distributed program, including test data, etc.: 1 799 469Distribution format: tar.gzProgramming language: MathematicaComputer: Any computer running unix and MathematicaOperating system: UnixSupplementary material: The Development guideClassification: 4.14, 5, 16.1, 20Nature of problem: Import data from output files of scientific computing packages, such as Gaussian, into Mathematica for symbolic calculation and production of publication quality tables and plots.Solution method: Provision of mnemonic top-down parsing procedures, functional programming.Running time: The complete extraction of data from a small basis density functional calculation on the water molecule, and from a larger basis density functional calculation on the zinc hydrate ion, that ran to 33 iterations, took 1 second and 23 seconds, respectively, on a Dell Poweredge 1750.     

Qprop: A Schrödinger-solver for intense laser-atom interaction

The Qprop package is presented. Qprop has been developed to study laser-atom interaction in the nonperturbative regime where nonlinear phenomena such as above-threshold ionization, high order harmonic generation, and dynamic stabilization are known to occur. In the nonrelativistic regime and within the single active electron approximation, these phenomena can be studied with Qprop in the most rigorous way by solving the time-dependent Schrödinger equation in three spatial dimensions. Because Qprop is optimized for the study of quantum systems that are spherically symmetric in their initial, unperturbed configuration, all wavefunctions are expanded in spherical harmonics. Time-propagation of the wavefunctions is performed using a split-operator approach. Photoelectron spectra are calculated employing a window-operator technique. Besides the solution of the time-dependent Schrödinger equation in single active electron approximation, Qprop allows to study many-electron systems via the solution of the time-dependent Kohn-Sham equations.

Program summary

Program title:QPROPCatalogue number:ADXBProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXBProgram obtainable from:CPC Program Library, Queen's University of Belfast, N. IrelandComputer on which program has been tested:PC Pentium IV, AthlonOperating system:LinuxProgram language used:C++Memory required to execute with typical data:Memory requirements depend on the number of propagated orbitals and on the size of the orbitals. For instance, time-propagation of a hydrogenic wavefunction in the perturbative regime requires about 64 KB RAM (4 radial orbitals with 1000 grid points). Propagation in the strongly nonperturbative regime providing energy spectra up to high energies may need 60 radial orbitals, each with 30000 grid points, i.e. about 30 MB. Examples are given in the article.No. of bits in a word:Real and complex valued numbers of double precision are usedNo. of lines in distributed program, including test data, etc.:69 995No. of bytes in distributed program, including test data, etc.: 2 927 567Peripheral used:Disk for input-output, terminal for interaction with the userCPU time required to execute test data:Execution time depends on the size of the propagated orbitals and the number of time-stepsDistribution format:tar.gzNature of the physical problem:Atoms put into the strong field of modern lasers display a wealth of novel phenomena that are not accessible to conventional perturbation theory where the external field is considered small as compared to inneratomic forces. Hence, the full ab initio solution of the time-dependent Schrödinger equation is desirable but in full dimensionality only feasible for no more than two (active) electrons. If many-electron effects come into play or effective ground state potentials are needed, (time-dependent) density functional theory may be employed. Qprop aims at providing tools for (i) the time-propagation of the wavefunction according to the time-dependent Schrödinger equation, (ii) the time-propagation of Kohn-Sham orbitals according to the time-dependent Kohn-Sham equations, and (iii) the energy-analysis of the final one-electron wavefunction (or the Kohn-Sham orbitals).Method of solution:An expansion of the wavefunction in spherical harmonics leads to a coupled set of equations for the radial wavefunctions. These radial wavefunctions are propagated using a split-operator technique and the Crank-Nicolson approximation for the short-time propagator. The initial ground state is obtained via imaginary time-propagation for spherically symmetric (but otherwise arbitrary) effective potentials. Excited states can be obtained through the combination of imaginary time-propagation and orthogonalization. For the Kohn-Sham scheme a multipole expansion of the effective potential is employed. Wavefunctions can be analyzed using the window-operator technique, facilitating the calculation of electron spectra, either angular-resolved or integratedRestrictions onto the complexity of the problem:The coupling of the atom to the external field is treated in dipole approximation. The time-dependent Schrödinger solver is restricted to the treatment of a single active electron. As concerns the time-dependent density functional mode of Qprop, the Hartree-potential (accounting for the classical electron-electron repulsion) is expanded up to the quadrupole. Only the monopole term of the Krieger-Li-Iafrate exchange potential is currently implemented. As in any nontrivial optimization problem, convergence to the optimal many-electron state (i.e. the ground state) is not automatically guaranteedExternal routines/libraries used:The program uses the well established libraries blas, lapack, and f2c     

A Mathematica program for the two-step twelfth-order method with multi-derivative for the numerical solution of a one-dimensional Schrödinger equation

In this paper, we present the detailed Mathematica symbolic derivation and the program which is used to integrate a one-dimensional Schrödinger equation by a new two-step numerical method. We add the fourth- and sixth-order derivatives to raise the precision of the traditional Numerov's method from fourth order to twelfth order, and to expand the interval of periodicity from (0,6) to the one of (0,9.7954) and (9.94792,55.6062). In the program we use an efficient algorithm to calculate the first-order derivative and avoid unnecessarily repeated calculation resulting from the multi-derivatives. We use the well-known Woods-Saxon's potential to test our method. The numerical test shows that the new method is not only superior to the previous lower order ones in accuracy, but also in the efficiency. This program is specially applied to the problem where a high accuracy or a larger step size is required.

Program summary

Title of program: ShdEq.nbCatalogue number: ADTTProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADTTProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputer for which the program is designed and others on which it has been tested: The program has been designed for the microcomputer and been tested on the microcomputer.Computers: IBM PCOperating systems under which the program has been tested: Windows XPProgramming language used: Mathematica 4.2Memory required to execute with typical data: 51 712 bytesNo. of bytes in distributed program, including test data, etc.: 45 381No. of lines in distributed program, including test data, etc.: 7311Distribution format: tar gzip fileCPC Program Library subprograms used: noNature of physical problem: Numerical integration of one-dimensional or radial Schrödinger equation to find the eigenvalues for a bound states and phase shift for a continuum state.Method of solution: Using a two-step method twelfth-order method to integrate a Schrödinger equation numerically from both two ends and the connecting conditions at the matching point, an eigenvalue for a bound state or a resonant state with a given phase shift can be found.Restrictions on the complexity of the problem: The analytic form of the potential function and its high-order derivatives must be known.Typical running time: Less than one second.Unusual features of the program: Take advantage of the high-order derivatives of the potential function and efficient algorithm, the program can provide all the numerical solution of a given Schrödinger equation, either a bound or a resonant state, with a very high precision and within a very short CPU time. The program can apply to a very broad range of problems because the method has a very large interval of periodicity.References: [1] T.E. Simos, Proc. Roy. Soc. London A 441 (1993) 283.[2] Z. Wang, Y. Dai, An eighth-order two-step formula for the numerical integration of the one-dimensional Schrödinger equation, Numer. Math. J. Chinese Univ. 12 (2003) 146.[3] Z. Wang, Y. Dai, An twelfth-order four-step formula for the numerical integration of the one-dimensional Schrödinger equation, Internat. J. Modern Phys. C 14 (2003) 1087.     

Optimized evaluation of a large sum of functions using a three-grid approach

The efficient evaluation of numerous values of functions that vary rapidly only in a small part of the region of interest is presented. An optimized algorithm using a sequence of grids with increasing resolution is developed. The algorithm does not make any assumptions about special properties of the function to be evaluated, e.g., symmetry. An additional speed-up is obtained by exploiting the asymptotic behaviour of the functions to be summed. Two applications from high resolution atmospheric radiative transfer modelling in the infrared and microwave are presented. In a third example asymmetric Rautian line shapes important for high resolution molecular spectroscopy are considered. Computational gains by more than two orders of magnitude with relative errors less than 10⁻³ have been achieved.