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.     

RADCAP: A potential model tool for direct capture reactions

A computer program is presented aiming at the calculation of bound and continuum states, reduced transition probabilities, phase-shifts, photo-disintegration cross sections, radiative capture cross sections, and astrophysical S-factors, for a two-body nuclear system. The code is based on a potential model of a Woods-Saxon, a Gaussian, or a M3Y, type. It can be used to calculate nuclear reaction rates in numerous astrophysical scenarios.

Program summary

Title of program: RADCAP (RADiative CApture)Catalogue identifier:ADSHProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADSHProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandComputers: The code has been created on an IBM-PC, but also runs on UNIX machinesOperating systems: WINDOWS or UNIXProgram language used: Fortran-77Memory required to execute with typical data: 8 Mbytes of RAM memory and 2 MB of hard disk spaceNo. of bits in a word: 32 or 64Memory required for test run with typical data: 2 MBNo. of bytes in distributed program: 376 817No. of lines in distributed program, including test data, etc.: 3054Distribution format: tar gzip fileKeywords: Potential model, photodissociation, radiative capture, astrophysical S-factorsNature of physical problem: The program calculates bound and continuum wavefunctions, phase-shifts and resonance widths, astrophysical S-factors, and other quantities of interest for direct capture reactions.Method of solution: Solves the radial Schrödinger equation for bound and for continuum states. First the eigenenergy is estimated by using the WKB method. Then, a Numerov integration is used outwardly and inwardly and a matching at the nuclear surface is done to obtain the energy and the bound state wavefunction with good accuracy. The continuum states are obtained by a Runge-Kutta integration, matching the Coulomb wavefunctions at large distances outside the range of the nuclear potential.Typical running time: Almost all the CPU time is consumed by the solution of the radial Schrödinger equation. It is about 1 min on a 1 GHz Intel P4-processor machine for a Woods-Saxon potential.     

A basis-set based Fortran program to solve the Gross-Pitaevskii equation for dilute Bose gases in harmonic and anharmonic traps

Inhomogeneous boson systems, such as the dilute gases of integral spin atoms in low-temperature magnetic traps, are believed to be well described by the Gross-Pitaevskii equation (GPE). GPE is a nonlinear Schrödinger equation which describes the order parameter of such systems at the mean field level. In the present work, we describe a Fortran 90 computer program developed by us, which solves the GPE using a basis set expansion technique. In this technique, the condensate wave function (order parameter) is expanded in terms of the solutions of the simple-harmonic oscillator (SHO) characterizing the atomic trap. Additionally, the same approach is also used to solve the problems in which the trap is weakly anharmonic, and the anharmonic potential can be expressed as a polynomial in the position operators x, y, and z. The resulting eigenvalue problem is solved iteratively using either the self-consistent-field (SCF) approach, or the imaginary time steepest-descent (SD) approach. Iterations can be initiated using either the simple-harmonic-oscillator ground state solution, or the Thomas-Fermi (TF) solution. It is found that for condensates containing up to a few hundred atoms, both approaches lead to rapid convergence. However, in the strong interaction limit of condensates containing thousands of atoms, it is the SD approach coupled with the TF starting orbitals, which leads to quick convergence. Our results for harmonic traps are also compared with those published by other authors using different numerical approaches, and excellent agreement is obtained. GPE is also solved for a few anharmonic potentials, and the influence of anharmonicity on the condensate is discussed. Additionally, the notion of Shannon entropy for the condensate wave function is defined and studied as a function of the number of particles in the trap. It is demonstrated numerically that the entropy increases with the particle number in a monotonic way.

Program summary

Title of program:bose.xCatalogue identifier:ADWZ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADWZ_v1_0Program obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandDistribution format:tar.gzComputers:PC's/Linux, Sun Ultra 10/Solaris, HP Alpha/Tru64, IBM/AIXProgramming language used:mostly Fortran 90Number of bytes in distributed program, including test data, etc.:27 313Number of lines in distributed program, including test data, etc.:28 015Card punching code:ASCIINature of physical problem:It is widely believed that the static properties of dilute Bose condensates, as obtained in atomic traps, can be described to a fairly good accuracy by the time-independent Gross-Pitaevskii equation. This program presents an efficient approach of solving this equation.Method of solution:The solutions of the Gross-Pitaevskii equation corresponding to the condensates in atomic traps are expanded as linear combinations of simple-harmonic oscillator eigenfunctions. Thus, the Gross-Pitaevskii equation which is a second-order nonlinear differential equation, is transformed into a matrix eigenvalue problem. Thereby, its solutions are obtained in a self-consistent manner, using methods of computational linear algebra.Unusual features of the program:None     

An object-oriented C++ implementation of Davidson method for finding a few selected extreme eigenpairs of a large, sparse, real, symmetric matrix

A C++ class named Davidson is presented for determining a few eigenpairs with lowest or alternatively highest values of a large, real, symmetric matrix. The algorithm described by Stathopoulos and Fischer is used. The exception mechanism is involved to report the errors. The class is written in ANSI C++, so it is fully portable. In addition a console program as well as a program with graphical user interface for Microsoft Windows is attached, which allow one to calculate the lowest eigenstates of time-independent Schrödinger equation for a given binding potential in one, two or three spatial dimensions. The package contains the classes providing often used potential functions (model atom potential, Coulomb potential, square well potential and Kramers-Henneberger well potential) as well as a possibility to use any potential stored in a file (then any dimensionality of the problem is allowed).The described code is the subject of M.Sc. thesis of T.D. prepared under the supervision of J.M.

Program summary

Program title: DavidsonCatalogue identifier: ADZM_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADZM_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.: 3 037 055No. of bytes in distributed program, including test data, etc.: 20 002 609Distribution format: tar.gzProgramming language: C++Computer: AllOperating system: AnyRAM: User's parameters dependentWord size: 32 and 64 bitsSupplementary material: Test results for the 2D and 3D cases is availableClassification: 4, 4.8Nature of problem: Finding a few extreme eigenpairs of a real, symmetric, sparse matrix. Examples in quantum optics (interaction of matter with a laser field).Solution method: Davidson algorithmRunning time: The test example included in the distribution package (1D matrix) takes approximately 30 minutes to run. 2D matrix calculations can take hours and 3D, days, to run.     

The two-electron atomic systems. S-states

A simple Mathematica program for computing the S-state energies and wave functions of two-electron (helium-like) atoms (ions) is presented. The well-known method of projecting the Schrödinger equation onto the finite subspace of basis functions was applied. The basis functions are composed of the exponentials combined with integer powers of the simplest perimetric coordinates. No special subroutines were used, only built-in objects supported by Mathematica. The accuracy of results and computation time depend on the basis size. The precise energy values of 7-8 significant figures along with the corresponding wave functions can be computed on a single processor within a few minutes. The resultant wave functions have a simple analytical form consisting of elementary functions, that enables one to calculate the expectation values of arbitrary physical operators without any difficulties.

Program summary

Program title: TwoElAtom-SCatalogue identifier: AEFK_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEFK_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.: 10 185No. of bytes in distributed program, including test data, etc.: 495 164Distribution format: tar.gzProgramming language: Mathematica 6.0; 7.0Computer: Any PCOperating system: Any which supports Mathematica; tested under Microsoft Windows XP and Linux SUSE 11.0RAM:?¹⁰⁹ bytesClassification: 2.1, 2.2, 2.7, 2.9Nature of problem: The Schrödinger equation for atoms (ions) with more than one electron has not been solved analytically. Approximate methods must be applied in order to obtain the wave functions or other physical attributes from quantum mechanical calculations.Solution method: The S-wave function is expanded into a triple basis set in three perimetric coordinates. Method of projecting the two-electron Schrödinger equation (for atoms/ions) onto a subspace of the basis functions enables one to obtain the set of homogeneous linear equations F.C=0 for the coefficients C of the above expansion. The roots of equation det(F)=0 yield the bound energies.Restrictions: First, the too large length of expansion (basis size) takes the too large computation time giving no perceptible improvement in accuracy. Second, the order of polynomial Ω (input parameter) in the wave function expansion enables one to calculate the excited nS-states up to n=Ω+1 inclusive.Additional comments: The CPC Program Library includes “A program to calculate the eigenfunctions of the random phase approximation for two electron systems” (AAJD). It should be emphasized that this fortran code realizes a very rough approximation describing only the averaged electron density of the two electron systems. It does not characterize the properties of the individual electrons and has a number of input parameters including the Roothaan orbitals.Running time: ∼10 minutes (depends on basis size and computer speed)     

P-stability, Trigonometric-fitting and the numerical solution of the radial Schrödinger equation

In this paper we will study the importance of the properties of P-stability and Trigonometric-fitting for the numerical integration of the one-dimensional Schrödinger equation. This will be done via the error analysis and the application of the studied methods to the numerical solution of the radial Schrödinger equation.     

MOMDIS: a Glauber model computer code for knockout reactions

A computer program is described to calculate momentum distributions in stripping and diffraction dissociation reactions. A Glauber model is used with the scattering wavefunctions calculated in the eikonal approximation. The program is appropriate for knockout reactions at intermediate energy collisions (30 MeV?Elab/nucleon?2000 MeV). It is particularly useful for reactions involving unstable nuclear beams, or exotic nuclei (e.g., neutron-rich nuclei), and studies of single-particle occupancy probabilities (spectroscopic factors) and other related physical observables. Such studies are an essential part of the scientific program of radioactive beam facilities, as in for instance the proposed RIA (Rare Isotope Accelerator) facility in the US.

Program summary

Title of program: MOMDIS (MOMentum DIStributions)Catalogue identifier:ADXZ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXZ_v1_0Computers: The code has been created on an IBM-PC, but also runs on UNIX or LINUX machinesOperating systems: WINDOWS or UNIXProgram language used: Fortran-77Memory required to execute with typical data: 16 Mbytes of RAM memory and 2 MB of hard disk spaceNo. of lines in distributed program, including test data, etc.: 6255No. of bytes in distributed program, including test data, etc.: 63 568Distribution format: tar.gzNature of physical problem: The program calculates bound wavefunctions, eikonal S-matrices, total cross-sections and momentum distributions of interest in nuclear knockout reactions at intermediate energies.Method of solution: Solves the radial Schrödinger equation for bound states. A Numerov integration is used outwardly and inwardly and a matching at the nuclear surface is done to obtain the energy and the bound state wavefunction with good accuracy. The S-matrices are obtained using eikonal wavefunctions and the “t-ρρ” method to obtain the eikonal phase-shifts. The momentum distributions are obtained by means of a Gaussian expansion of integrands. Main integrals are performed with the Simpson's method.Typical running time: Almost all CPU time is consumed by calculations of integrals, specially for transverse momentum distributions which involves multiple integral loops. It takes up to 30 min on a 2 GHz Intel P4-processor machine.     

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 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.     

ALTDSE: An Arnoldi–Lanczos program to solve the time-dependent Schrödinger equation

We describe a general ab initio and non-perturbative method to solve the time-dependent Schrödinger equation (TDSE) for the interaction of a strong attosecond laser pulse with a general atom. While the field-free Hamiltonian and the dipole matrices may be generated using an arbitrary primitive basis, they are assumed to have been transformed to the eigenbasis of the problem before the solution of the TDSE is propagated in time using the Arnoldi–Lanczos method. Probabilities for survival of the ground state, excitation, and single ionization can be extracted from the propagated wavefunction.

Program summary

Program title: ALTDSECatalogue identifier: AEDM_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEDM_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.: 2154No. of bytes in distributed program, including test data, etc.: 30 827Distribution format: tar.gzProgramming language: Fortran 95. [A Fortran 2003 call to “flush” is used to simplify monitoring the output file during execution. If this function is not available, these statements should be commented out.].Computer: Shared-memory machinesOperating system: Linux, OpenMPHas the code been vectorized or parallelized?: YesRAM: Several Gb, depending on matrix size and number of processorsSupplementary material: To facilitate the execution of the program, Hamiltonian field-free and dipole matrix files are provided.Classification: 2.5External routines: LAPACK, BLASNature of problem: We describe a computer program for a general ab initio and non-perturbative method to solve the time-dependent Schrödinger equation (TDSE) for the interaction of a strong attosecond laser pulse with a general atom [1,2]. The probabilities for survival of the initial state, excitation of discrete states, and single ionization due to multi-photon processes can be obtained.Solution method: The solution of the TDSE is propagated in time using the Arnoldi–Lanczos method. The field-free Hamiltonian and the dipole matrices, originally generated in an arbitrary basis (e.g., the flexible B-spline R-matrix (BSR) method with non-orthogonal orbitals [3]), must be provided in the eigenbasis of the problem as input.Restrictions: The present program is restricted to a ¹S^e initial state and linearly polarized light. This is the most common situation experimentally, but a generalization is straightforward.Running time: Several hours, depending on the number of threads used.References: [1] X. Guan, O. Zatsarinny, K. Bartschat, B.I. Schneider, J. Feist, C.J. Noble, Phys. Rev. A 76 (2007) 053411. [2] X. Guan, C.J. Noble, O. Zatsarinny, K. Bartschat, B.I. Schneider, Phys. Rev. A 78 (2008) 053402. [3] O. Zatsarinny, Comput. Phys. Comm. 174 (2006) 273.     

Computation of the eigenvalues of the Schrödinger equation by exponentially-fitted Runge-Kutta-Nyström methods

In this work we consider exponentially fitted and trigonometrically fitted Runge-Kutta-Nyström methods. These methods integrate exactly differential systems whose solutions can be expressed as linear combinations of the set of functions exp(wx), exp(−wx), or sin(wx), cos(wx), w∈ℜ. We modify existing RKN methods of fifth and sixth order. We apply these methods to the computation of the eigenvalues of the Schrödinger equation with different potentials as the harmonic oscillator, the doubly anharmonic oscillator and the exponential potential.     

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.     

Application of the Feynman-Kac path integral method in finding excited states of quantum systems

Group theory considerations and properties of a continuous path are used to define a failure tree procedure for finding eigenvalues of the Schrödinger equation using stochastic methods. The procedure is used to calculate the lowest excited state eigenvalues of eigenfunctions possessing anti-symmetric nodal regions in configuration space using the Feynman-Kac path integral method. Within this method the solution of the imaginary time Schrödinger equation is approximated by random walk simulations on a discrete grid constrained only by symmetry considerations of the Hamiltonian. The required symmetry constraints on random walk simulations are associated with a given irreducible representation and are found by identifying the eigenvalues for the irreducible representation corresponding to symmetric or antisymmetric eigenfunctions for each group operator. The method provides exact eigenvalues of excited states in the limit of infinitesimal step size and infinite time. The numerical method is applied to compute the eigenvalues of the lowest excited states of the hydrogenic atom that transform as Γ² and Γ⁴ irreducible representations. Numerical results are compared with exact analytical results.     

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.     

The solution of stationary ODE problems in quantum mechanics by Magnus methods with stepsize control

In solid state physics the solution of the Dirac and Schrödinger equation by operator splitting methods leads to differential equations with oscillating solutions for the radial direction. For standard time integrators like Runge-Kutta or multistep methods the stepsize is restricted approximately by the length of the period. In contrast the recently developed Magnus methods allow stepsizes that are substantially larger than one period. They are based on a Lie group approach and incorporate exponential functions and matrix commutators. A stepsize control is implemented and tested. As numerical examples eigenvalue problems for the radial Schrödinger equation and the radial Dirac equation are solved. Further, phase shifts for scattering solutions for hydrogen atoms and copper are computed.     

TiReX: Replica-exchange molecular dynamics using Tinker

We present a driver program for performing replica-exchange molecular dynamics simulations with the Tinker package. Parallelization is based on the Message Passing Interface, with every replica assigned to a separate process. The algorithm is not communication intensive, which makes the program suitable for running even on loosely coupled cluster systems. Particular attention is paid to the practical aspects of analyzing the program output.

Program summary

Program title: TiReXCatalogue identifier: AEEK_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEEK_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.: 43 385No. of bytes in distributed program, including test data, etc.: 502 262Distribution format: tar.gzProgramming language: Fortran 90/95Computer: Most UNIX machinesOperating system: LinuxHas the code been vectorized or parallelized?: parallelized with MPIClassification: 16.13External routines: TINKER version 4.2 or 5.0, built as a libraryNature of problem: Replica-exchange molecular dynamics.Solution method: Each replica is assigned to a separate process; temperatures are swapped between replicas at regular time intervals.Running time: The sample run may take up to a few minutes.     

g_permute: Permutation-reduced phase space density compaction

Biomolecular processes are governed by free energy changes and thus depend on a fine-tuned interplay between entropy and enthalpy. To calculate accurate values for entropies from simulations is particularly challenging for the solvation shell of proteins, which contributes crucially to the total entropy of solvated proteins, due to the diffusive motion of the solvent molecules. Accordingly, for each frame of a Molecular dynamics (MD) trajectory, our software relabels the solvent molecules, such that the resulting configuration space volume is reduced by a factor of N! with N being the number of solvent molecules. The combinatorial explosion of a naive implementation is here overcome by transforming the task into a linear assignment problem, for which algorithms with complexity O(N³) exist. We have shown in previous research that the solvent entropy can be estimated from such a compacted trajectory by established entropy estimation methods. In this paper, we describe the software implementation which also allows applications beyond entropy estimation, such as the permutation of lipids in membrane bilayers.

Program summary

Program title: g_permuteCatalogue identifier: AECJ_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AECJ_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: GPLNo. of lines in distributed program, including test data, etc.: 45 173No. of bytes in distributed program, including test data, etc.: 2 730 678Distribution format: tar.gzProgramming language: CComputer: PC-compatible running LinuxOperating system: LinuxRAM: Dependent on the number of solvent molecules, min 12 582 912 bytesClassification: 3, 4.8, 4.9External routines: liblap (included); From GROMACS-3.3.1: libgmx (not included)Nature of problem: Estimating the entropy of solvent molecules from a molecular dynamics simulation trajectory cannot be performed on ordinary trajectories.Solution method: Compacting the configuration space of molecules by exploiting their permutation symmetry. Applies to trajectories either compatible to those obtained with the GROMACS simulation package [1] or multi-model pdb (Protein Data Bank) files.Restrictions: In rare cases the time to find a solution for the linear assignment problem can be very long.Running time: Dependent on trajectory length and number of molecules to be permuted.References:[1] D. van der Spoel, et al., J. Comput. Chem. 26 (2005) 1701.